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1.0  INTRODUCTION 


In  support  of  the  Navy's  continuing  program  to  develop  high  speed 
towed-body  systems,  this  report  describes  the  extension  of  an  earlier 
analysis  capabilities  programming  effort.  It  is  the  results  from  this 
most  recent  work  which  we  are  reporting  here. 

To  provide  some  background  for  the  present  tasks,  a  brief  review  of 
our  previous  efforts  is  given  below. 

The  initial  work  was  undertaken  to  provide  an  analysis  and  numerical 
demonstration  capability  for  the  study  of  streamlined  fairings.  These 
fairings  are  to  be  used  in  high  speed  underwater  towing  situations.  The 
construction  of  the  fairings  is  such  that  streamwise  distortions  are  able 
to  occur;  and,  that  a  likely  consequence  from  them  would  be  some  "unde¬ 
sirable  behavior"  from  the  towline.  What  was  not  anticipated,  but  what 
has  actually  occured,  was  an  uncontrolled  "kiting"  of  these  tow  cables  — 
leading  to  unacceptable  towing  performance. 

Since  the  physics  of  this  problem  was  not  immediately  apparent,  and 
the  remedial  measures  taken  to  correct  the  adverse  happenings  did  not 
resolve  the  situation;  it  was  proposed  that  some  numerical  analysis  tools 
be  developed  to  aid  in  the  study  and  understanding  of  these  situations. 

A  first  step  in  this  effort  was  the  development  of  the  computational 
algorithm  which  we  have  Tabled  HYDROSAP.  This  program  seeks  to  describe 
the  physical  deformations  produced  on  a  typical  streamlined,  flexible, 
two-dimensional  fairing  section  as  a  consequence  of  the  hydrodynamic 
loads  impressed  on  it.  The  loadings  which  these  airfoil-type  sections 
are  subjected  to  are  those  typical  to  any  streamlined  (two-dimensional) 
shape  operating  at  an  angle  of  attack  in  a  steady,  viscous  fluid  flow. 
Because  these  fairings  are  flexible  —  their  after-sections  are  con¬ 
structed  from  a  rubber-like  material  —  the  chordwise  distortions  cause 
the  flow  field  to  be  altered;  which,  in  turn,  produces  other  chordwise 
distortions;  which,  in  turn,  changes  the  flow  pattern,  •••  etc. 
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As  a  consequence  of  these  actions  it  is  apparent  that  the  mathematical 
modelling  of  this  phenomenon  requires  an  iterative  method  of  solution.  A 
solution  which  (first)  determines  the  hydrodynamic  "loads"  on  the  profile 
section;  then  ascertains  the  distortions  due  to  this  loading;  subse¬ 
quently,  redefining  the  loads  on  the  (deformed)  profile;  and  so  on. 
(Incidentally,  this  methodology  is  depicted,  schematically,  on  Figure  2.1 
in  the  next  section  of  this  report.) 

Being  cognizant  of  the  fact  the  HYDROSAP  treats  only  the  two- 
dimensional  problems;  that  it  will  merely  "point"  to  the  probable  causes 
of  kiting;  then,  a  next  logical  step  (in  the  developmental  process)  was 
that  of  modelling  the  three-dimensional  problem  and  seeking  a  more  likely 
and  realistic  answer  to  the  real-world  situation. 

The  three-dimensional,  computational  program  development  is  basi¬ 
cally  what  is  being  reported  on  here.  However,  it  should  be  noted  that 
we  have  only  begun  the  task  of  developing  a  three-dimensional,  numerical 
analysis  tool.  In  this  report  we  will  describe  the  current  status  of 
these  efforts;  outline  what  has  been  accomplished;  discuss  the  findings 
and  results  from  this  work;  and  describe  those  things  which  we  believe 
are  needed  to  accomplish  the  goals  for  this  overall  task.  These  goals, 
incidentally,  will  only  be  realized  when  a  full,  three-dimensional  anal¬ 
ysis  program  has  been  developed;  and  after  it  has  been  checked  and  veri¬ 
fied  against  experimental  evidence.  Then,  and  only  then,  can  this  work 
be  said  to  be  "complete". 

In  this  second  stage  of  the  development  effort  two  primary  computa¬ 
tional  algorithms  have  been  developed.  One,  given  the  acronym  SHCENT,  is 
used  to  describe  the  mechanical-structural  characteristics  of  the  cable's 
cross-sections.  Its  principal  outputs  are  the  shear-center  location  and 
torsional  rigidity;  both  of  which  depend  on  the  shape  and  the  material 
description  of  the  cross-sections.  In  addition,  certain  other  parameters 
(elastic  and  geometric  properties)  needed  in  the  structural  analysis  are 
also  computed  here. 


With  the  loading  information  provided  from  HYDROSAP,  and  the  cross- 
sectional  mechanical  properties  obtained  from  SHCENT,  the  second  computa¬ 
tional  program,  TOWLINE,  is  used  to  ascertain  the  three-dimensional 
dynamic  response  of  the  cable.  This  algorithm  was  developed  using  a 
finite  element  method  technique.  The  cable's  elements  were  build  around 
a  standard  beam  model,  one  which  is  assuned  to  suffer  large  displacements 
but  have  only  small  strains. 

The  TOWLINE  program,  in  its  present  state  of  development,  is  des¬ 
cribed  in  this  document.  It  should  be  recognized  that  the  existing  ver¬ 
sion  of  TOWLINE  is  not  ready  to  be  used  on  a  production  basis,  even 
though  it  has  been  used  in  an  exploratory  fashion  on  several  "model" 
cable  problems  as  we  describe  in  Section  3.0,  below.  The  original  intent 
of  this  developmental  task  was  to  provide  the  existing  computational 
framework;  to  enhance  it  during  subsequent  efforts;  and,  to  perfect  the 
model  through  a  verification  with  experimental  data  and  evidence. 

In  the  following  sections  of  this  document  we  will  (first)  outline 
our  previous  developments  (HYDROSAP);  then,  we  will  describe,  discuss  and 
comment  on  the  achievements  and  results  obtained  in  the  current  task 
effort. 
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2.0  TWO-DIMENSIONAL  FLOW-STRUCTURE  INTERACTION 


2.1  Overview  of  the  Two-Dimensional  Problem 


A  study  of  the  nonlinear  interactions  between  hydrodynamic  and 
structural  effects,  on  a  flexible  two-dimensional  streamlined  cross- 
section  when  exposed  to  a  steady,  viscous  fluid  flow,  was  a  first  step  in 
the  analysis  of  faired  towing  cable  systems.  The  resulting  cross- 
section's  deformed  shape  alters  the  cable's  local  hydrodynamic  character¬ 
istics  and,  consequently,  affects  its  global  behavior.  This  build-up  in 
interaction  of  structural  and  hydrodynamic  effects  is  pictured  schematic¬ 
ally  in  Figure  2.1. 

In  our  iterative  approach  to  this  problem  we  determine  a  pressure 
distribution  corresponding  to  the  base  (or  design)  profile  shape.  Thi: 
streamlined  section  is  assumed  to  be  at  a  fixed  angle  of  attack,  a  , 
relative  to  the  free  stream  velocity  vector,  V  .  This  (base)  pressure 
distribution  is  then  applied  to  a  suitably  described  finite  element  model 
of  the  fairing  to  determine  a  structurally  "deformed”  profile.  This 
altered  shape,  in  turn,  develops  an  altered  pressure  distribution  -- 
which  in  turn  deforms  the  profile  —  which  in  turn  etc.  This  pro¬ 
cess  continues,  in  an  iterative  fashion,  until  the  incremental  change  in 
structural  deformation  is  "sufficiently  small"  to  describe  the  converged 
(deformed)  profile  and  its  attendant  pressure  distribution.  Typically  a 
converged  solution  is  achieved  after  some  five  to  ten  iterations. 

The  HYDROSAP  computer  program  is  a  mechanization  of  the  process 
described  above.  HYDROSAP  is  based  on  various  elements  of  two  existing 
computer  programs;  the  NASA  Multicomponent  Two-dimensional  Viscous 
Airfoil  Program  (referred  to  as  HYDRO  in  the  sequel)  and  the  University 
of  California's  Non-linear  Structural  Analysis  Program  (NONSAP).  The 
various  program  elements,  plus  the  many  components  developed  for  our  pro¬ 
grams’s  mechanizations,  were  integrated  into  a  single  software  system. 
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{To  a  converged  solution) 

Figure  2.1 

A  Schematic  Depicting  the  Iterative  Operations  of  HYDROSAP  -- 
A  Hydrodynamic  Loading,  Structural  Response  Software  System 
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In  sections  2.2.1  and,  respectively,  2.2.2  below,  we  will  briefly 
describe  these  solution  algorithms.  A  more  detailed  discussion, 
including  a  description  of  program  input  and  operating  requirements,  can 
be  found  in  the  HYDROSAP  documentation  [l].1 

That  report  also  describes,  in  full  detail,  the  several  features  of 
the  HYDROSAP  system.  In  particular,  a  description  of  the  BTS-developed 
executive,  DRIVER,  and  other  features  which  would  be  of  interest,  primar¬ 
ily,  to  the  software  user,  appear  therein. 

2.2  Solution  Algorithms 

2.2.1  The  Fluid  Flow  Field 

Given  a  pointwise  shape  description  of  a  two-dimensional  profile, 
the  HYDRO  routines  compute  the  velocity  destribution,  over  the  profile's 
surface,  by  a  superposition  of  velocity  fields.  These  are  the  primary 
velocity  distributions  due  to  thickness  and  to  profile  camber.  This  is  a 
potential  flow  solution,  one  which  must  be  modified  to  account  for  vis¬ 
cous  (boundary  layer)  effects  by  one  of  two  methods:  (1)  a  modifying  of 
the  profile's  geometry;  or,  (2)  introducing  a  distribution  of  source  sin¬ 
gularities  to  account  for  the  boundary  layer  thickness  distribution.  The 
potential  flow  field  (thickness  and  camber  effects)  is  then  re-computed, 
taking  into  account  the  boundary  layer  displacement  thickness  over  the 
profile.  This  process  reoccurs,  iteratively,  until  changes  in  the  veloc¬ 
ity  field  fall  below  some  pre-set  tolerance.  Of  course,  if  the  viscous 
boundary  layer  is  neglected,  the  computations  are  not  iteratively 
sequenced  and  the  corresponding  run  time  is  considerably  shortened. 

Calculations  for  the  velocity  fields,  due  to  thickness  and  camber, 
make  use  of  Oellers'  Vortex  Distribution  Method  [5],  using  the  scheme 
which  is  outlined  in  [2]  and  [3].  Incorporated  into  these  computations 
is  a  procedure  to  satisfy  the  Kutta  condition  for  the  profile.  The  pro¬ 
cedure  represents  a  modified  Kutta  condition  wherein  the  upper  and  lower 


Numerals  in  brackets  [  ]  denote  the  reference  nimbers  in  Section  5.0. 
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surface  vortex  strengths  are  matched,  in  magnitude,  but  assigned  opposing 
signs  at  the  foil's  downstream  terminus. 

If  the  viscous  boundary  layer  is  to  be  simulated  by  a  distribution 
of  source  singularities  then  subroutine  SOURCE  is  called.  Otherwise  the 
boundary  layer's  influence  is  provided  for  by  modifying  the  velocity 
field  for  camber  effects.  For  those  cases  where  the  foil's  trailing  edge 
thickness  is  not  zero,  the  velocity  field  due  to  thickness  is  altered 
through  a  closure  of  the  geometry  one  chord  length  downstream  of  the  true 
trailing  edge. 

The  present  version  of  HYDRO  does  not  compute  flow  separation, 
although  recent  work  in  this  area  could  be  incorporated  [6].  Also,  the 
drag  increment  due  to  the  wake  behind  the  foil  is  based  solely  on  empir¬ 
ical  results  from  tests  on  truncated  profiles. 

The  laminar  boundary  layer  computations  in  HYDRO  are  based  on 
Goradia's  adaptation  of  the  work  by  Cohen  and  Reshotko  [7],  and  on 
Schlichting  [8].  In  addition  to  computing  the  laminar  boundary  layer 
thickness,  the  program  determines  the  laminar  friction  coefficients  and 
predicts  boundary  layer  transition  points.  Truckenbrodt1  s  integral 
method  [9]  is  used  to  compute  parameters  describing  the  turbulent  bound¬ 
ary  layer. 

Finally,  the  usual  force  and  moment  coefficients  for  the  profile  are 
computed  from  the  pressure  and  shear  force  distributions.  The  moment 
coefficients  —  at  both  the  quarter-chord  point  and  at  the  leading  edge 
—  are  determined.  In  addition  to  the  "usual"  drag  coefficient,  the  drag 
coefficient  based  on  the  approximate  theory  of  Squire  and  Young  [10]  is 
calculated  and  displayed  in  the  output. 
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2.2.2  A  Cross-Section's  Mech-anical  Response 


Given  a  description  of  the  profile's  topology,  the  mechanical  prop¬ 
erties  of  its  constituent  materials,  and  the  pressure  loading  (from 
HYDRO),  NONSAP  is  used  to  compute  the  profile's  mechanical  response. 

(The  finite  element  program  NONSAP  (NON! inear  Structural  Analysis 
Program)  was  originally  developed  in  1974  [12]).  Even  though  NONSAP  may 
be  used  to  analyze  the  static  or  dynamic  response  of  one-,  two-,  or 
three-dimensional  structures,  including  the  effects  of  geometric  and/or 
material  nonlinearities,  we  restrict  ourselves  to  its  two-dimensional 
capabilities. 

In  this  section  we  give  a  brief  overview  of  NONSAP's  capabilities 
and  its  operation.  However,  for  a  more  complete  and  detailed  description 
of  the  NONSAP  theory  and  operation,  the  reader  is  referred  to  the  orig¬ 
inal  program  documentation  [11,  12]  and  to  the  HYDROSAP  documentation 
[1]. 

2. 2. 2.1  Program  Description 

The  NONSAP  program  was  developed  to  solve  static  or  dynamic,  linear 
or  nonlinear  structural  problems.  The  allowable  nonlinearities  may  be 
due  to  either  material  behavior  (elastic,  hyperelastic,  and  hypoelastic 
material  models  are  available)  or  to  the  effects  of  large  displacements 
and  large  strains.  The  equations  of  motion  from  continuum  mechanics, gov¬ 
erning  the  behavior  of  the  structure  under  study,  are  discretized,  in 
space  variables,  through  the  use  of  isoparametric  finite  elements.  The 
equations  are  discretized  in  the  time  variable  through  either  a  Wilson- 
Theta  or  a  Newmark-Beta  time  integration  schane.  The  resulting  nonlinear 
equations  are  solved  using  an  incremental  method  which  corresponds  to  a 
modified  Newton  iteration. 

The  structural  systems  to  be  analyzed  may  be  modeled  using  a  variety 
of  finite  elanent  types.  These  are: 
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•  one-dimensional  truss  elements, 

•  two-dimensional  plane  stress  or  plane  strain  elements, 

•  two-dimensional  axi- symmetric  shell  elements,  and 

•  three-dimensional  solid  or  thick  shell  elements. 

Since  nonlinearities  may  be  due  to  nonlinear  material  behavior  as 
well  as  to  large  displacements,  or  large  strains,  there  are  several 
material  descriptions  available  in  the  program.  These  include: 

•  an  isotropic  linear  elastic  material  model, 

•  an  orthotropic  linear  elastic  material  model, 

•  the  Mooney-Rivl in  model, 

•  an  elastic-plastic  material  with  von  Mises  or  Drucker-Prager 
yield  conditions, 

•  a  variable  tangent  model,  and 

•  a  curve  description  model. 

Currently,  the  HYDROSAP  program  makes  use  only  of  the  one-dimensional 
truss  finite  element,  the  two-dimensional  plane  strain  element,  and  the 
isotropic  linear  elastic  material  model.  All  nonlinearities  are,  there¬ 
fore,  geometric.  Moreover,  HYDROSAP  uses  NONSAP  routines  in  the  non¬ 
linear  static  solution  mode  with  the  "updated  Lagrangian"  solution  proce¬ 
dure.  As  a  consequence  of  the  constraints,  the  discussions  below  will  be 
restricted  to  those  NONSAP  procedures  and  options  which  are  currently 
employed  in  HYDROSAP. 

2. 2. 2. 2  Incremental  Equilibrium  Equations  for  Structural  Systems 

Using  the  well-known  and  well  documented  [12,  13]  method  of  finite 
elements,  the  incremental  nodal  point  equilibrium  equations  for  an 
assemblage  of  finite  elements  take  the  form. 


M  u(t+At)  +  c  u( t+At)  +  K(t)  Au  =  R(t+At)  -  F(t) ,  (2.1) 


where 


M  =  the  constant  mass  matrix, 

C  =  the  constant  damping  matrix, 

K(t)  =  the  tangent  stiffness  matrix,  described  at  time  t  , 

R(t)  =  the  external  applied  load  vector,  described  at  time  t  , 

F(t)  =  the  nodal  point  force  vector  (of  stress-equivalent  internal 
reactions),  for  time  t  , 

u(t)  =  the  vector  of  nodal  point  displacements,  at  time  t  , 

Au  =  the  vector  of  nodal  point  displacement  increments,  from  time 

t  to  time  t+At:  thus,  u(t+At)  -  u(t)  , 
t  =  time,  and 

At  =  time  increment. 

Of  course,  in  a  static  analysis  M=C=0  ,  while  time  is  regarded  only  as  a 
loading  parameter  used  for  applying  the  external  load  in  "quasistatic" 
load  increments.  Also,  in  nonlinear  analysis,  the  solution  of  (2.1) 
yields  only  approximate  solution  increments,  Au  .  In  order  to  improve 
the  solution  accuracy  and  to  prevent  the  buildup  of  solution  instabil¬ 
ities,  equilibriim  iterations  may  be  used  at  selected  time  steps. 

In  the  nonlinear  static  case  an  equation  system  of  the  form 

K(t)(Au).+1  =  R(t+At)  -  F.  (i=0,l,2,  •••)  (2.2) 

is  solved  by  iteration.  Here  FQ  =  F(t)  and  F.  is  the  vector  of 
stress-equivalent  internal  reactions  corresponding  to  the  displacement 
vector  u.j (t+At)  =  u.  ^(t+At)  +  (Au)^  .  The  solution  of  (2.2),  by 
iteration,  is  clearly  a  modified  Newton  method  for  (2.1). 
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2. 2. 2. 3  NONSAP  Solution  Process 


The  complete  NONSAP  solution  process  consists  of  three  distinct 
stages:  input,  assembly,  and  step-by-step  solution  (See  Figure  2.2). 

Input  Stage 

In  this  stage,  control  and  nodal  point  input  are  read  and/or  gener¬ 
ated  by  the  program;  and  equation  numbers  for  the  active  degrees  of 
freedom  at  each  nodal  point  are  established.  In  addition,  externally 
applied  load  vectors  are  established  and  stored  on  tape.  Finally,  when 
the  elenent  data  are  ready  and/or  generated,  element  connectivity  arrays 
are  established,  and  all  element  data  are  stored  on  tape. 

Assembly  Stage 


Here  the  linear  and  effective  linear  structural  stiffness  matrices 
are  assembled  and  stored  on  tape,  as  are  the  mass  and  damping  matrices 
(for  a  dynamic  analysis).  Matrices  are  stored  in  compacted  form  as  one¬ 
dimensional  arrays.  Only  the  elements  below  the  "skyline"  of  a  matrix 
are  processed,  thereby  reducing  storage,  input/output,  and  computational 
costs. 

Solution  Stage 

In  nonlinear  static  analysis  (the  solution  mode  for  which  NONSAP 
routines  are  used  by  HYDROSAP)  damping  and  mass  effects  are  neglected. 
Time  serves  as  a  loading  parameter  only,  and  the  input  time  step  deter¬ 
mines  the  load  increnent  to  be  applied  to  the  structure  during  each  solu¬ 
tion  step.  The  linear  stiffness  matrix,  which  was  previously  assembled 
and  stored  on  tape,  is  updated  at  each  solution  step.  This  is  accom¬ 
plished  by  updating  the  stiffness  matrices,  of  the  nonlinear  elements,  to 
form  a  current  tangent  stiffness  matrix.  The  load  step  interval  at  which 
this  update  occurs  is  one  of  the  NONSAP  control  inputs. 
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Figure  2.2 

fiONSAP  Flowchart  Program  Operation 
in  Nonlinear  Static  Solution  Mode 


Moreover,  within  each  load  step,  solution  accuracy  may  be  improved 
by  equilibrium  iterations.  The  interval  of  load  steps  at  which  an  equi¬ 
librium  iteration  is  to  be  allowed  is  also  a  NONSAP  control  input. 

A  solution  to  the  incremental  linear  equations  at  each  load  step  is 
performed  by  the  linear  equation  solver,  COLSOL,  which  performs  Gauss 
eliminations  on  the  positive  definite  symmetrical  system  of  equations. 
This  algorithm  takes  advantage  of  the  banded  and  sparse  structure  of  the 
stiffness  equations  by  not  processing  elements  outside  a  matrix  skyline, 
since  they  remain  zero  throughout  the  computation.  The  algorithm 
consists  of  an  LDL^  decomposition  of  the  tangent  stiffness  matrix  into 
lower  triangular,  L  ,  and  diagonal,  D  ,  matrices,  followed  by  a  reduc¬ 
tion  and  back  substitution  of  the  load  vector. 

2. 2. 2. 4  NONSAP  Element  Library 

In  this  section  we  will  briefly  describe  only  those  NONSAP  element 
types  which  are  used  by  HYDROSAP. 

The  One-Dimensional  Truss  Element 


A  one-dimensional  truss  element,  which 
may  be  located  in  three-dimensional  space, 
is  used  in  HYDROSAP  to  model  the  thin  skin 
which  is  on  the  exterior  of  some  cable  fair¬ 
ings.  This  element  has  two  nodes;  it  is 
assumed  to  have  a  constant  cross-sectional 
area  (which  corresponds  to  the  skin  thickness 
in  the  HYDROSAP  plane  strain  analysis)  and  to 
undergo  small  strains  and  large  displace¬ 
ments. 
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Figure  2.3 
A  One-Dimensional 
Truss  Element 


The  Two-Dimensional  Plane  Stress/Plane  Strain  Element 


An  isoparametric  plane  strain  element, 
having  a  variable  number  of  nodes,  is  used  to 
model  the  cross-section  of  the  cable  fairing. 
This  element  may  have  from  three  to  eight 
nodes;  any  of  the  nodes  numbered  from  five 
through  eight  on  the  sketch  may  be  omitted.  A 
three  node  element  is  obtained  by  setting  the 
coordinates  of  two  nodes  for  a  four  node  quad¬ 
rilateral  equal  to  one  another.  In  a  plane 
strain  analysis  the  out-of-plane  element  thick¬ 
ness  is  assumed  to  be  unity. 


x 


Figure  2.4 
A  Two-Dimensional 
Plane  Strain  Isopara 
metric  Quadrilateral 
Element 


2.3  Sample  Computations 

In  order  to  illustrate  the  capabilities  of  the  HYDROSAP  software 
system,  in  this  section  we  will  present  the  results  of  some  typical  com¬ 
putations.  In  section  2.3.1  will  we  analyze  the  behavior  of  an  NACA  0020 
section.  In  section  2.3.2  we  will  discuss  an  NACA  65-018  section;  and  in 
section  2.3.3  we  will  explore  the  sensitivity  of  the  computations  to  a 
refinement  of  the  finite  element  mesh. 

2.3.1  The  NACA  0020  Profile  Section 


A  hypothetical,  three-inch  chord  faired  cable,  constructed  as  an 
NACA  0020  section,  is  "towed"  through  water  at  44  fps  (34.5  knots).  It 
is  assumed  to  be  at  an  angle  of  attack  of  8°.  In  the  following  discus¬ 
sion  we  designate  parameters  of  the  computer  model  by  capitalization. 
These  variables  are  defined  in  Appendix  A. 
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Structurally,  the  towline  section  is  modelled  as  having  a  rigid 
CORE  (ECORE  =  lxlO6  psi)  extending  aft  from  the  leading  edge  for  3/4" 
(Y1C0RE  =  0.0,  Y2C0RE  =  0.75).  The  afterbody,  or  TAIL,  is  assumed  to  be 

4 

composed  of  a  relatively  flexible  rubber-like  material  (ETAIL  =  1x10 
psi).  The  entire  fairing  section  is  enclosed  by  a  very  flexible  skin 
(ESKIN  =  2x10^  psi),  which  is  0.05  inches  thick  ( TSKIN  =  0.05).  The 
finite  element  mesh  used  to  analyze  the  structural  response  is  generated 
automatically  by  the  software  (MESHC  =  11,  MESHT  =  3).  This  mesh  con¬ 
sists  of  48  two-dimensional  elements  having  57  nodes;  the  outer  SKIN  con¬ 
necting  the  surface  nodes  is  modelled  by  24  one-dimensional  elements. 

For  purposes  of  the  hydrodynamic  computations,  the  fairing's  NACA 
0020  surface  profile  is  modelled  by  35  surface  points  distributed  auto¬ 
matically  (HYDROSAP  controls:  IPANEL  =2,  and  N  =  -2  •••  for  F0IL03)  on 
the  upper  and  lower  fairing  surfaces  by  the  software.  Viscous  effects 
are  included  in  the  computations  in  each  of  the  HYDROSAP  "outer  iter¬ 
ations"  ( ITRSWT  =  1). 

The  effect  of  the  hydrodynamic  loading  on  this  flexible  symmetric 
fairing  is  to  induce  negative  camber.  Figure  2.5  summarizes  the  types  of 
deformations  suffered  by  the  fairing  cross-section.  Evidently  the  mater¬ 
ial  above  the  fairing  chord  is  in  compression  while  that  below  is  in  ten¬ 
sion.  The  equilibrium  trailing  edge  transverse  deflection  is 
approximately  2.3%  of  the  chord. 

The  negative  camber  acts  to  reduce  the  magnitude  of  the  lift  and 
(nose)  moment  coefficients  significantly.  In  Figure  2.6  these  hydro- 
dynamic  coefficients  are  plotted  for  each  successive  pass  through  the 
HYDROSAP  system.  It  is  evident  from  the  graph  that  after  ten  iterations 
the  hydrodynamic  loading  and  structural  response  are  in  equilibrium.  The 
fairing's  deformation  has  reduced  the  lift  coefficient  by  24%,  from  = 
0.8996  for  the  undeformed  fairing,  to  =  0.6794  for  the  deformed 
fairing.  The  nose-down  moment  coefficient,  taken  about  the  profile's 
leading  edge,  has  been  reduced  by  29%.  Figure  2.7  shows  the  effect  that 
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Figure  2.5 

Transverse  Deflection  of  Centerline 
for  a  Flexible  Cable  Fairing 
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Figure  2.7 

Pressure  Distributions  Over 
Undeformed  and  Deformed  Fa i 


the  modest  chordwise  deflection  of  the  fairing  has  on  the  overall  pres¬ 
sure  distribution.  Pressure  coefficients  are  reduced  in  magnitude  by 
about  10%  to  20%  even  in  the  region  of  the  rigid  leading  edge.  The 
trailing  edge  deflection  evidently  makes  its  presence  felt  globally  over 
the  entire  fairing's  surface. 

2.3.2  The  NACA  65-018  Profile  Section 


We  explore  the  degree  of  deformation  for  a  flexible  NACA  65-018 
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fairing  made  entirely  of  a  material  with  Young's  modules  (ETAIL)  =  1x10 
psi.  No  distinguishable  SKIN  or  CORE  are  present  here.  The  three-inch 
section  is  towed  at  44  fps  (34.5  kts)  through  sea-water. 

In  Figure  2.8  we  plot  the  viscous  and  potential  pressure  distribu¬ 
tions  over  the  undeformed  section  at  0°  angle  of  attack.  Since  the  sec¬ 
tion  is  symmetric  the  upper  and  lower  surface  pressure  distributions  are 
identical,  and  the  profile  suffers  no  deformations.  Figures  2.9  and  2.10 
show  pressure  distributions  for  the  deformed  and  undeformed  fairing  at  2° 
and  8°  angles  of  attack,  respectively.  Figure  2.11  summarizes  the  angle 
of  attack  dependence  of  the  sectional  lift,  drag,  and  nose  moment  coef¬ 
ficients. 

Since  the  NACA  65-018  fairing  material  is  ten  times  stiffer  than 
that  of  our  NACA  0020  example,  the  development  of  negative  camber  is 
somewhat  less  pronounced.  Nevertheless,  the  gradual  separation  of  the 
deformed  and  undeformed  lift  and  moment  curves  in  Figure  2.11  underscores 
the  need  for  careful  attention  to  be  given  to  the  local  orientation  and 
behavior  of  the  cable  sections  in  a  three-dimensional  analysis.  We  note 
that  the  drag  coefficients  are  not  so  sensitive  to  the  pressure  induced 
chordwise  deformation  of  the  fairing. 
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Figure  2.8 

Pressure  Distribution  for  Potential  and  Viscous  Fluid  Flows 
Over  a  Symmetrical  (NACA  65-018)  Profile  at  0°  Angle  of  Attack 
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and  Deformed  Fairing  at  2°  Angle  of  Attack 
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2.3.3  An  Accuracy  Study 

Although  certain  theoretical  error  estimates  are  available,  to  meas¬ 
ure  the  accuracy  of  a  particular  finite  element  solution,  in  practice  it 
is  customary  to  determine  the  sensitivity  of  the  numerical  solution  to 
the  fineness  of  the  finite  element  mesh  by  experiment.  In  this  regard  we 
compute  the  numerical  solution  for  a  sequence  of  finite  element  meshes; 
each  succeeding  mesh  being  finer  than  its  predecessor. 

For  this  study  we  use  the  NACA  0020  fairing  of  section  2.3.1,  but 
tow  it  at  4°  angle  of  attack,  at  22  fps  (13  kts),  in  seawater.  We  use 
the  number  of  nodal  points  for  the  mesh  as  a  measure  of  mesh  "fineness”. 
For  this  analysis  the  solutions  were  computed  using  meshes  with  7,  27, 

57,  79,  and  90  node  points.  Table  2.1  summarizes  the  conditions  of  the 
tow  and  the  nominal  hydrodynamic  coefficients  for  the  rigid  profile. 

Table  2.2  gives  the  variation  in  trailing  edge  deflection,  the  hydro- 
dynamic  coefficients,  and  computer  “costs"  with  the  fineness  of  the 
mesh.  Evidently,  increasing  the  number  of  nodes  from  79  to  90  offers  a 
negligible  increase  in  accuracy.  Indeed,  even  the  crudest  (and 
"stiffest")  finite  element  mesh  (having  only  17  nodes)  already  accounts 
for  25%  of  the  reduction  in  ,  observed  for  the  deformed  fairing 
modelled  with  90  nodes.  Obtaining  the  remaining  15%  reduction 
requires  an  investment  of  twice  as  much  computer  CPU  time. 

This  completes  the  brief  review  of  our  two-dimensional  analysis 
capabilities.  In  the  next  section  a  description  of  the  three-dimensional 
analysis  development  will  be  presented. 
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Table  2.1 


Nominal  Test  Values 


PROFILE 

NACA  0020  section 

CHORD  (length) 

3.0  (in.) 

CORE  (length) 

0.5  (in.) 

TAIL  MATERIAL  MODULUS 

104  (psi.) 

TEST  CONDITIONS 

V 

22.  (fps.) 

a 

+4.  (deg.) 

HYDRODYNAMIC  COEFFS. 

C£ 

0.4507 

Sno 

-0.1223 

Cd 

0.0086 

Table  2.2 

Sensitivity  Analysis  Values 


No.  of 
Nodes 

Trailing  Edge 
Deflection 

Coefficients1 

No.  of 
Iterations^ 

Run  Time3 
CPU  I/O 

(mins) 

17 

-0.0016  0.0084 

0.4248  -0.1141  0.0085 

3 

0.50  4.70 

27 

0.0086 

0.4243  -0.1139  0.0085 

3 

57 

0.0097 

0.4205  -0.1126  0.0084 

4 

0.90  5.62 

79 

0.0098 

0.4203  -0.1125  0.0084 

4 

1.07  6.06 

90 

0.0098 

0.4202  -0.1125  0.0084 

4 

1.15  6.09 

NOTTS:  ’Coefficients  determined  from  Pressure  Distributions; 
"Outer  Loop  Iterations,  to  a  coverged  solution; 

3Based  on  an  IBM  360/91  system. 
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3.0  THE  THREE-DIMENSIONAL  FLOW-STRUCTURE  INTERACTION 


3.1  Overview  of  the  Three-Dimensional  Problem 


In  contrast  to  our  focus  on  the  local  deformations  for  a  typical 
cable  cross-section,  in  the  two-dimensional  study,  the  object  of  the 
three-dimensional  study  is  to  determine  the  global  shape  of  a  tow  cable. 
In  this  case  the  cable  is  being  influenced  by  forces  from  the  towed  body 
in  addition  to  the  hydrodynamic  forces  developed  on  its  cross-sectional 
profiles.  For  this  study  the  solution  methodology  is  based  on  a  large 
deformation,  small  strain  theory  for  the  deflections  of  a  heterogeneous, 
isotropic  rod.  These  deflections  are  modelled  in  terms  of  the  three- 
dimensional  motion  of  a  reference  curve  —  the  locus  of  shear  centers  for 
the  cable's  cross-sections.  The  cross-section's  mechanical  properties, 
such  as  bending  stiffnesses,  torsional  rigidity,  modulus  weighted  cen¬ 
troids,  and  shear  center  locations  are  the  subject  of  a  separate  anal¬ 
ysis;  this  is  detailed  in  section  3.2.1,  below. 

Having  reduced  the  motion  of  the  towing  cable  to  that  of  the  locus 
of  shear  centers,  the  dynamic  equations  for  the  displacements  of  the  ref¬ 
erence  curve  are  discretized  using  the  finite  element  technique.  We  dis¬ 
cuss  this  process  and  its  computer  implementation  in  section  3.2.2. 

3.2  Solution  Algorithms 

3.2.1  Mechanical  Properties  for  a  Cable  Cross-Section 

In  order  to  model  a  towing  cable's  dynamical  behavior,  in  terms  of 
the  motion  of  its  reference  axis,  it  is  necessary  to  compute  for  a  typ¬ 
ical  cross-section  (or  cross-sections,  in  the  case  of  a  non-uniform 
cable)  certain  of  its  mechanical  properties. 
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We  wish  to  compute  the  following  properties  for  any  representative 
section: 

•  moduli  3  weighted  area,  "EA  =  /  E  dy  dz 

Q 


•  torsional  rigidity, 


GJ  =  /  G(y2+z2+yf ,  -zf.J  dy  dz  . 
a  y 

Here  the  function  f  is  the  St.  Venant  warping  function  for  the  sec¬ 
tion.  The  warping  function  describes  the  longitudinal  (out  of  plane) 
deformation  for  the  cross-secion,  considered  to  be  a  thin  prismatic  mem¬ 
ber  subject  to  pure  torsion. 


By  way  of  explanation,  consider  the 
prismatic  member  shown  in  Figure  3.2;  let  it 
be  subject  to  a  twist  per  unit  length  of 
along  the  axis  of  the  member.  The  0  dis¬ 
placement  field,  for  pure  torsion,  is 
described  as: 


x 


Figure  3.2 
A  General  Prismatic 
Member  Under  Twist,  0 

u  =  0  f(y,z) 

v  =  -0  x  z  (3.1) 


w  =  0  x  y  , 

where  (u,v,w)  are  displacements,  due  to  O  ,  in  the  (x,y,z)  directions. 
For  displacements  of  the  form  (3.1),  the  non-zero  strains  are: 


and 


yx  3x 


3v  +  iM  =  0  (1!  _z) 

3y  0  l3y  z) 


zx 


is  +  =  0  (”  +y) 

3x  3z  v3z 


(3.2) 


The  equations  of  equilibrium,  for  a  heterogeneous,  isotropic  material, 
therefore  reduce  to 


(3.3a) 


in  fl  ,  subject  to  the  boundary  condition 


(V'i').n  =  zny  -  ynz 


(3.3b) 


on  a  . 

In  equation  (3.3b),  n  =  (n^,  nz)  is  the  outward  normal  to  . 
Clearly,  (3.2)  and  (3.3)  determine  f  only  up  to  a  constant;  thus  we 
will  normalize  ¥  by  requiring  that 

/  G  V  dy  dz  =  0  .  (3.4) 

n 

The  computer  program  SHCENT  (SHear  CENTer)  is  used  to  compute  the 
warping  function  (i*)  for  the  airfoil-type  cross-sections  used  for 
faired  towlines  whose  towing  characteristics  are  to  be  studied.2 


2SHCENT  is  an  adaptation  of  a  code  developed  by  E.  M.  E.  Raumgartner 
[15]. 
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The  finite  element  method  is  used  to  solve  equations  (3.3).  For 
this,  the  cable's  cross-section  is  triangulated;  ¥  is  assumed  to  be 
linear  on  each  element  of  the  triangul ation ;  the  shear  modulus,  G  ,  is 
assumed  constant  over  each  element  of  the  triangulation;  and  a  linear 
equation  for  the  nodal  values  of  ’f  is  obtained  from  the  variational 
form  of  equations  (3.3).  We  outline  this  procedure  below. 

Let  4>  be  any  function,  defined  on  the  cross-section  ft  ,  which 
together  with  its  gradient  is  square  integrable  over  ft  .  Then  the  weak 
(or  variational)  form  of  equations  (3.3)  is 

-2)  +  G  (|f  +y)  If]  dy  dz  =  0  .  (3.5) 

To  obtain  equation  (3.5)  we  have  multiplied  (3.3a)  by  <j>  ,  integrated 
over  ft  ,  applied  Green's  Theorum,  and  finally  taken  advantage  of  (3.3b). 


Suppose,  for  the  moment,  that  ft  is  th:  ». 

triangle,  T  ,  as  in  Figure  3.3,  and  that  we  /  3\ 

consider  v  to  be  of  the  form  /  j  \ 

(ypzj  (y2’Z2) 


Figure  3.3 

Sketch  of  a  Triangular 
Element  with 
Coordinates  (y,z). 

i  (y.z)  =  2I  (fi  ^(y.z)  +  \  f2(y,z)  +  t3  f3(y,z)J  (3.6) 
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where 


^  (y.z)  s  ( y-Yj )  (Zj-zk)  -  (z-Zj)  (yj-yk)  ,  (3.7) 

for,  respectively, 

i=l,2, 3 
j=2,3,l 
k=3,l ,2 

and  where  A  is  the  area  of  T  .  Vie  take  <j>  in  analogous  form.  Then 
clearly  f  (y.  ,z^)  =  ^  ,  i=l,2,3,  and  f  is  linear  on  T  . 


Next  let  us  observe  that 


Similar  equations,  of  course,  hold  for  <P  •  Thus,  equation  (3.5) 
becomes 


L  2A  V23TV2z3i’ Vl2J  3*  _ t ,'ly23 ^ T2  y31*'3yl?'  3z}dy 


/  G  (^  -  y||  )  dy  dz  . 
fi  J 


(3.10) 


Since  equation  (3.10)  must  hold  for  arbitrary  <f>  ,  we  select,  alterna¬ 
tively,  4>  =  f. ,  i=l  ,2,3.  Thus,  we  have  from  equation  (3.10)  three 
equations  for  the  three  unknowns,  » ^2 * ^3  *  *n  ma^r1x  ^orri  Those  three 
equations  are 


K  *  =  F 


(3.11) 


where:  ^(Vj.^.f^T,  F  sfFj.Fg,  F2)T,  and  K  =  (k...)  for  i,j=l,2,3. 
The  matrix  K  is  called  the  "local  stiffness  matrix";  its  elements  are 
given  by 

kj-^  =  G(y22+Z23)/2A  , 


k12  =  k21  G^y23y31+Z23Z31^2A  ’ 


k13  _  kj^  =  G  (y^^y, 5+z9^z1 ?)/2A  , 


2ZJ  12  23  12y 


(3.12) 


k22  =  ^3 (y + z )/2A  , 


31  31  ■ 


and 


k23  "  k32  =  G^y31y12+Z31Z12^2A  * 


k33  =  G^y12+Z12^2A 
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The  vector  of  "generalized  loads"  is  given  by 


F1  =  *23+i  z23^  * 

F2  =  G(y  y31+z  z31)  ,  (3.13) 


and 


F3  =  G (y  y12+z  z12)  . 


Here  (y  ,  z)  are  the  coordinates  of  the  centroid  of  the  element,  T  . 

We  retain  the  shear  coefficient,  G  ,  in  both  equation  (3.12)  and  equa¬ 
tion  (3.13)  since  the  global  stiffness  matrix  and  global  generalized  load 
vector,  for  the  entire  domain  fi  ,  must  be  assembled  from  the  local 
values  for  each  element.  The  local  values  of  the  shear  constant  and  ele¬ 
ment  area  may  vary  from  element  to  element. 

The  global  analog  of  equation  (3.11);  i.e.,  the  fully  assembled 
matrix  problem,  is  a  linear  equation  for  the  nodal  values  of  the  "best" 
piecewise  linear  approximation  to  T  for  the  given  triangulation  of 
[12].  The  global  problem  is  solved  by  elimination,  taking  advantage, 
however,  of  the  banded,  symmetric  and  sparse  structure  of  the  global 
stiffness  matrix. 

3. 2. 1.1  Torsional  Rigidity 

The  torsional  rigidity,  GO  ,  of  the  cross-section  Q  is  defined  to 
be  the  moment  necessary  to  twist  a  rod  (of  cross-section  n  )  one  radian 
per  unit  of  length.  That  is,  if  the  twist  of  the  rod  is  0  ,  then  the 
necessary  twisting  moment,  M  ,  is 


Setting  0=1,  and  writing  M  in  terms  of  the  cross-sectional  shear 
stresses,  we  have 


1 


gj  =  /  (-zaVv+ya,*  )  dy  dz 

-  G  -z(|fz)  +  y(-|Ify)  dy  dz 

=  /  G  (y2+z2+yf,  -zY,  )  dy  dz  .  (3.15) 

a  L  y 

Several  remarks  concerning  equation  (3.15)  are  in  order. 

If  G  is  constant  over  .ft  ,  we  see  that 
_  2  2 

GJ  =  G*  /  (y  +z  +yY  -zY ,  )  dy  dz  .  (3.16) 

ft  y 

The  integral  on  the  right  hand  side  of  (3.16)  depends  only  on  the  shape 
of  ft  ,  and  is  called  the  torsion  constant  —  usually  denoted  by  J  — 
for  the  cross-section.  When  ft  is  a  circular  section,  f  =  0  over  ft 
then  J  reduces  to  this  section's  polar  moment  of  inertia. 

3.2. 1.2  Bending  Rigidities 


The  bending  rigidities,  EI^  ,  EIzz  ,  and  EIyz  are  computed  in 
SHCENT  by  calculating  the  bending  rigidities  for  each  triangular  element 
T  ,  about  its  own  centroid.  Using  the  transfer  of  axes  theorem,  these 
element  rigidities  are  referred  to  the  modulus  weighted  centroid  of  the 
section.  Bending  rigidities  with  respect  to  the  global  origin  --  the 
modulus-weighted  centroid  —  and  with  respect  to  the  shear  center  are 
among  the  program' s  output. 


3.2. 1.3  Modulus  Weighted  Area  and  Centroid 

As  with  the  bending  rigidities,  the  modulus  weighted  area  and 
weighted  centroid,  for  a  ,  are  (each)  computed  by  summing  over  all  the 
triangular  elements,  T  ,  into  which  ft  is  partitioned;  thus: 


TK  =  l  E,  AT  ; 
T  1  1 


(3.17) 


also 


y  = 


f  *T  EI  AI 

rk 


and 


z 


j  ZT  ET  AT 

TK 


(3.18a) 


(3.18b) 


3. 2. 1.4  Shear  Center 


The  shear  center  for  fl  is  defined  to  be  that  point,  in  the  plane  of 
n  ,  at  which  a  transverse  point  load  P  =  (P  ,Pz)  ,  on  a  rod  of  cross- 
section  G  ,  causes  no  twist  about  its  longitudinal  axis.  It  is  well- 
known  [16]  that  the  resulting  normal  stress,  ,  at  a  point  (y,z) 
in  a  section,  at  distance  L  from  the  loaded  section,  is 

°xx  =  E(y,z)(C2y+C2z)L  ,  (3.19a) 
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where 


C,  = 


P  El  +  P  El 

=  y  .yy _ 2, _ xl. 


1  '  e r 


yy 


EIzz  -  ^ 


(3.19b) 


P  El  +  P  El 
r  =  v  yz  z  zz 

l2  ‘  " 


El 


yy 


EI 


zz 


try' 


(3.19c) 


Using  the  technique  in  reference  [17],  it  is  possible  to  show  that  the 
coordinates,  (y  ,  z$)  of  the  shear  center  are: 


yc  =  Cl  /  E(y,z)  yf(y.z)  dy  dz 
5  1  n 

+  Cl  /  E(y,z)  Zf(y,z)  dy  dz  ,  (3.20a) 

L  n 


and 


z.  =  C"  /  E(y,z)  yf(y,z)  dy  dz 
s  i  n 

+  C"  /  E(y,z)  zf(y,z)  dy  dz  .  (3.20b) 

^  n 

Here  the  constants  Cl  are  obtained  from  equations  (3.19)  when  Py  =  0  , 
Pz  =  1  ;  and  ,  C'l  are  the  values  for  Py  =  1  and  =  0  .  In  the 
SHCENT  computer  program,  the  integrals  in  equations  (3.20)  are  carried 
out  piecewise  over  n  ,  making  use  of  the  constancy  of  E(y,z)  and  the 
linearity  of  f(y,z)  over  each  element  T  in  the  triangulation  of  fl  . 


A  short  user's  guide  for  the  SHCENT  program  is  included  in  Appendix 
B  of  this  report.  Included  in  this  user's  guide  is  a  description  of  the 
ptogram' s  input  deck,  external  file  requirements,  and  a  flow  chart. 

3.2.2  The  Towline  Response 


We  formulate  the  equations  of  motion  for  the  discretized  towline 
in  terms  of  the  incremental  displacements  of  the  structure  as  it  moves 
from  one  (quasi)  equilibrium  configuration  to  an  adjacent  deformed  con¬ 
figuration.  If  we  denote  by  B(t)  the  deformed  state  of  the  towline  at 
time  t  ,  then  we  are  interested  in  the  incremental  deformation  "AB"  , 
of  the  towline,  from  state  B(t)  to  state  B(t+At)  .  We  compute  the 
displacements  which  define  AB  with  respect  to  a  local  connected  coordi¬ 
nate  system  using  a  linearized  theory.  These  local  deformations  are  then 
referred  to  a  fixed  global  coordinate  system  by  means  of  local  transfor¬ 
mations  which  depend  on  B(t)  . 

Let  us  drop,  for  the  moment,  the  dependence  on  t  and  focus  on  the 
computation  of  the  incremental  displacement. 

As  shown  in  Figure  3.4,  the  local  incre¬ 
mented  displacement  field  for  a  prismatic 
member,  B  ,  with  cross-section  fi(x)  ,  may 
be  written  in  terms  of  the  incremental  dis¬ 
placements  of  a  reference  axis.  That  is,  the 
displacements  of  material  points  in  body  B 
can  be  approximated,  over  the  cross-section, 

by 


AB  . 


z 

C^T  * 


Figure  3.4 
A  General  Prismatic 
Member 


I 


u(x,y,z)  s  u0(x)  -  <J>(x)z  -  ^(x)y  +  O' (x)  i*(x,y,z)  ,  (3.21a) 

v(x,y,z)  5  vQ(x)  -zO(x)  ,  (3.21b) 


and 


w(x,y,z)  s  wQ(x)  +  yo(x)  , 


(3.21c) 


where  u0«v0»w0  are  displacements  at  the  reference  axis;  0  is  the 
angle  of  twist;  4>  =  dwQ/dx  ,  and  V  =  dvQ/dx  are  the  bending  about  the 
y  and  z  axes,  respectively.  Also,  here,  ¥  is  the  St.  Venant  warping 
function  for  the  cross-section;  it  is  obtained  by  solving  the  St.  Venant 
torsion  problem  for  fl(x)  ,  as  discussed  in  section  3.2.1,  above.  In  the 
sequel,  we  assume  that  does  not  vary  along  the  length  of  the  member. 

The  strains  corresponding  to  the  displacement  field  (3.21),  and 
subject  to  the  assumption  ¥  =  ¥(y,z)  ,  are 


e  =  4^  =  u'  -  z<f>*  -  yT'  +  Y0"  , 
xx  ax  o  y  J  ’ 


3V  n 

eyy  “  sy  “  0  ’ 


9W  0 

ezz  =  3Z  =  °  * 


e  =  4^  +  4^  =  V  -  Z0'  -  'F  +  Y,  , 
xy  ax  ay  o  y  ’ 


_aw  _av  _ 
eyz  "ay  3z  ’ 


(3.22a) 

(3.22b) 

(3.22c) 

(3.22d) 

(3.22e) 


and 
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(3.2 2f) 


zx 


3u  _3w 
3z  3x 


-  4>  +  0'  f,  +  w'  +  yO' 

Z  0 


Here  the  prfmes  signify  differentiation  with  respect  to  x  .  We  relate 
the  non-zero  stresses  to  the  strains  through  the  constitutive  relations: 

axx  =  (X+2G)exx  ,  (3.23a) 

o  =  G  e  ,  (3.23b) 

xy  xy 

and 

a  =  G  e  .  (3.23c) 

zx  zx 

Here  X  =  vE/(l+v)(l-2v)  is  the  material's  Lame'  constant,  v  is  the 
Poisson  ratio,  E  is  the  Young's  modulus,  and  G  =  E/2(l+v)  is  the 
shear  modulus.  We  will  adopt  the  notation,  E  =  X  +  2G,  henceforth. 

The  internal  strain  energy  increment,  V  ,  of  body  B  ,  subject  to 
the  incremental  strain  field  (3.22),  and  incremental  stress  field  (3.23), 
is 


v  •  \  4  *  g  4 +  G  4 } d,<  d* dz  •  <3-24> 

Substituting  expressions  (3.22)  into  (3.24),  and  integrating  over 
the  cross-section,  gives  the  strain  energy  in  terms  of  the  six  functions 
u(x) ,  v(x) ,  w(x) ,  e(x),  <|>(x) ,  and  f(x)  ,  and  the  several  mechanical 
properties  of  the  cable's  cross-section: 

•  modulus  weighted  area,  lA  =  /  E  dy  dz  ; 

si 
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modulus  weighted  centroid 


y  =  ==  /  y  E  dy  dz  , 
EA  a 


and 


z  =  ~  z  /  E  dy  dz  ; 
EA  n 


•  bending  stiffnesses, 

El  =  /  z2  E  dy  dz  , 
y  a 

El  =  /  y2  E  dy  dz  , 
z  a 


and 


tlV7  =  /  yz  f  dy  dz 

O 


and , 

•  torsional  rigidity, 


GJ  =  /  GjyV+yT.^zt.yj  dy  dz  . 


We  have  discussed  the  computation  of  these  section  properties  at 
length  in  section  3.2.1.  There  are  several  other  cross-sectional 
properties  which  are  of  interest;  we  will  denote  these,  generically, 
using  the  notation: 


I(«)  =  /  <*(y,z)  dy  dz  .  (3.25) 

n 


Thus,  for  example,  we  could  write 


EA  =  1(E)  , 


or 


EIZ  =  I(Ey2)  . 

Using  this  notation,  and  integrating  over  n  in  equation  (3.24),  we 
obtain: 


•u 


EA  ((u^  -  2_zu;w"  -  2yu^)  +  2I(E>)u^ 


+  2  EIyz  u^  -  2I(Ezf)  6TW;  -  2I(Eyf)o;v^ 

+  EIy  (Wq)2  +  EIZ  (v^)2  +  I(Ef2)(0^)2 

+  GJ  (O'  )2}  dx  .  (3.26) 


-42- 


In  order  to  discretize  the  equations  of  motion  for  body  B  ,  we 
shall  assume  that  it  has  been  sub-divided  into  a  finite  number  of  para¬ 
metric  elements  (which  we  will  again  denote  by  B  ).  Within  each  of 
these  elements  the  incremental  displacements  (of  the  reference  curve)  can 
be  modelled  as  polynomials  in  the  connected  coordinate  system.  The  coef¬ 
ficients  of  the  polynomials  will  be  called  the  generalized  displacements, 
associated  with  that  element.  We  will  model  the  transverse  displace¬ 
ments,  v  and  w  ,  as  third  order  polynomials  in  the  axial  coordinate; 
but  the  twist,  0  ,  and  the  elongation,  u  ,  will  be  linear.  Note  that 
for  convenience  we  have  dropped  the  zero  subscripts  on  u,v,w.  Thus,  we 


assume: 


u(x)  =  Ujttj (x)  +  u2a  2 (x) 

v(x)  =  VjB^x)  +  v|Tj(x)  +  v262(x)  +  v^y2(x) 


(3.27a) 


(3.27b) 


w(x)  =  w^U)  +  w^y^x)  +  w2B2(x)  +  w^U)  .  (3.27c) 


O(x)  =  ©jOjfx)  +  02a2(x)  ; 


(3.27d) 


where,  with  L  being  the  length  of  the  element,  the  basis  (or  shape) 
functions  are 


al(x)  =  1  -  p 


a2(x)  =  J-  , 


0j(x)  =  l-3(£)2+  2(f)3 


(3.28a) 


(3.28b) 


(3.28c) 


Yjtx)  =  L|f  2{^)2  +  (£)3},  (3.28d) 

32(x)  =  l-3(l-£)2  +  2(l-£  )3  ,  (3.28e) 


and 


Y2(x)  =  -  L  2+  (1-1  )3J.  (3.28f ) 

(Recall  that  x  is  the  longitudinal  coordinate  for  the  element).  We 
note  that  since  v'  (x)  =  f(x)  and  w'(x)  =  <t>(x)  ,  we  may  write 


v 


i 


1-1,2 


and 


1-1,2  . 


By  virtue  of  the  definitions  (3.28)  —  for 
the  basic  functions  —  we  see  that  the 
generalized  displacements  (u.,  v^ ,  w.. ,  0. , 
fj  ,  for  i=l,2)  are  the  values  of  the  shape 
functions  (3.27)  at  the  node  points  (points  1 
and  2  in  Figure  3.5)  which  define  the  loca¬ 
tion  of  the  ends  of  the  "beam".  As  a  conse¬ 
quence  of  (3.27)  and  (3.28)  we  see  that 
0"=0  ,  and  so  the  internal  energy  becomes,  in 
terms  of  the  generalized  displacement, 


Figure  3.5 
A  Prismatic  Element 
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( U«_  U.  )  U/j-U.  r\ 

V  =  EA—^ - z  EA(-^)Z(^-^)-  y 


u2"ul, 


+  EIyz|(vrv2)C^|(wrw2)  +^(f1+4'2)] 

+  (i}>j+<}>2)C-^2(Wj-w2)  + -^[y2(  i(i^+2()>2)  +  ^  (2<t>^+4>2)]| 


£i 


+  ~/  P(wl"w2)2  +p(wrw2)(V‘*2)  +  I<*l+*2+*l*2)j 


+  ~r  {f§(vrv2>2  *  ^f(,r,2)(V’,2)  * 


_  (02-Oj)‘ 

♦ 1 


(3.29) 


Let  us  denote  by  |>l  =  (u^,  u2,  ^1’  ^1*  ^2*  ^2’  ^1*  ^2*  ^2^ 

the  vector  of  generalized  incremental  displacements.  We  see,  then,  that 
the  internal  energy  can  be  expressed  as  a  quadratic  form  in  the 
incremental  generalized  displacements;  or. 


v  -IrV- 


(3.30) 


where  K  is  the  local  element  stiffness  matrix  given  in  Table  3.1. 


Table  3.1 

A  12  x  12  Beam  Stiffness  Matrix 


Corresponding  to  there  is  a  vector,  P  ,  which  describes  the  net 
generalized  force  increments.  Vector  P  is  composed  of  the  net  force  on 
body  B  ;  i.e.,  P  is  determined  by  the  difference  between  the  exter¬ 
nally  applied  loads  for  the  configuration  B+AB  and  the  internal 
resisting  loads,  corresponding  to  the  deformed  configuration  B  .  We 
model  the  net  forces  on  B  as  loads  distributed  linearly  over  the  length 
of  the  element: 


Pr  '  Pr1al(x)  +  Pr2a2(x)  • 


and 


mr  =  +  mr  °^(x)  * 


(3.31a) 


(3.31b) 


Here  the  subscript  r  represents  x,  y,  or  z  ;  and  P^  ,  m  are  the 
force  and  moment  per  unit  length  acting  on  B  at  the  reference  axis. 
The  work  of  the  external  applied  loads  is  given  by 


W  =  /  ( Px  u+  Py v+  P2 why  0+iy  (j>+mz  v )  dx  , 


(3.32) 


where  u,  v,  w,  0,  <f>,  and  ¥  are  given  by  (3.27).  In  terms  of  the 
components  P  and  fj.  ,  W  is  given  by 


“  4  px,“l  *  Vpx,u2*px2ul>  *  7px?u2 


*  kj*l  4  (n,x®2™x20P  *  K292 


+ 


17L  31  L 

2rpyivl  +  20pyiv2  "  30 py 


¥ 

1  2 
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+ 


+ 


+ 


+ 

The  external  work. 


between  the  vector, 
of  generalized  loads 
by: 


lffPy2vl  +  Iff  py2fl  +  lrpy2v2 


17L  3L  Vr 

W\*l  +  20PZlW2  '  30Pz^2 


lfepZoWl  +  30Pz  *1  +  lrpz,w2 


(m  «nv  )(w2-w.)  (m  -m  )(VVL 
yl  y2  1  1  ,  y2  yl  *  1 

- -  + - r5 - 


(in 


,  -Hn  )(v2-v  )  (m  -m 

Z1  Z2  L  1  ,  z2__l____ 


12 


(3.33) 


I  ,  can  therefore  be  written  as  the  scalar  product, 


W  *  |tTP 


(3.34) 


|x  ,  of  generalized  displacements,  and  a  vector,  P  , 
,  Now,  for  clarity,  P  2  (Pj , • **,P12)T  is  defined 


L  , 


L  , 
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(3.35) 


-P..  L 


12 


~W 


+ - 12 


It  remains  to  compute  the  "kinetic  energy",  T  ,  of  the  displacement 
increment  from  state  R  to  state  B+AB  ;  i.e., 

T  -  \l  p(x,y,z) j(u)2  +  (v)2  +  (w)2jdx  dy  dz  .  ( 3 . 36 ) 3 

The  substitution  of  equations  (3.26)  into  (3.33),  with  the 
application  of (3.27)  and  (3.28),  results  in  an  equation  of  the  form: 

T  =  \  |iTMcji  (3.37) 

where  the  consistent  mass  matrix,  Mc  ,  is  a  result  of  the  integration  of 
the  material  density,  p  ,  against  the  various  products  of  shape  functions 
arising  as  a  consequence  of  (3.36).  We  note  that  Mc  is  a  full 
matrix.  In  computer  implementations  of  the  method  it  is  customary  to 
replace  M  by  a  diagonal  lumped  mass  matrix,  M  ,  in  the  interests  of 
speeding  the  numerical  integration  of  the  equations  of  motion.  The 
diagonal  elements  of  this  matrix,  M  ,  are  given  by: 


"l.l 

pAL 

=  2  » 

^2,2 

pAL 

=  2  * 

pLI 

^3,3 

X 

"  2  * 

3The  superscript  dot  denotes  a  time  derivative;  (*)  =  (*) 
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"4,4 

pLIx 

2  • 

M5,5 

pAL 

‘  2  ’ 

M6,6 

pAL3 

24  + 

PUZ 

2 

M7,7 

pAL 

2  * 

M8,8 

pAL3 

24 

PLIZ 

2  ’ 

^9,9 

pAL 

2  ’ 

.,3  pH 

Mio,io =  +  ~r  •  ^3*38^ 

The  kinetic  energy,  T  ,  potential  energy,  V  ,  and  external  work,  W  , 
for  the  entire  structure,  are  assembled  from  the  local  values  by  means  of 
the  local  to  global  element  transformation  matrices,  .  Using  the 
subscript  B  for  elemental  quantities,  and  the  subscript  G  for  global 
quantities,  we  may  write,  formally: 

T  =  Bz  2Fb  eb  ^  ^Pb’TPg  Vg  ’ 

V  =  J  2  f*B  EB  KB  EB  ^B  5  2  ^G  KG  P’G’ 
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and 


W  =  g  EB  EP  PB  H  Pg  P  '  (3.39) 

In  the  sequel  we  suppress  the  subscript  G  ;  thus,  all  generalized 
vectors  and  matrices  will  be  global  quantities  unless  otherwise  noted. 

Let  us  define  the  Lagrangian,  L  ,  for  this  deformation  increment  to 
be 

L  *  T-V+W  .  (3.40) 

Applying  Hamilton's  Principle,  we  require  that 

t+At 

J()i)  =/  (T-V+W)  ds  (3.41) 
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1 


be  extremized  by  the  solution  vector,  p  .  Setting  3J/8(a=  0  gives 
Langrange's  equations: 


i_  ik  Ik 

dt  d\i  3K 


(3.42) 


However,  since  we  know  that 


and 


d_ 

dt 


then  equation  (3.42)  implies  that 


Mjl+  kjx  =  P  .  (3.43) 

We  integrate  equation  (3.42)  using  the  implicit  Newmark  beta  method 
[11]  as  follows:  Introducing  the  notation, 

p(t)  =  piQ  ,  (=0  ,  from  the  definition  of  p.  )  , 


|U(t+At)  =px  , 


p(t)  =p0  , 


|i(t+At)  =  ^  , 

f^(t)  =p.g  , 
and 

ji(t+At)  =ji1  , 

We  express  the  displacement  increment,  and  velocity,  at  the  end  of  the 
time  step  by 

h  =  Ho  +  Ho4t  +  (i-s)poAt2  +  <Vt?  « 

and 

•  Mo+Mi 

Mi  =Po  +c¥1  At  *  (3.44) 

We  seek  ^  ^  and  |ij,  given  fxQ,  jiQ  and  |i0(=0).  Premulti¬ 

plying  the  first  of  equations  (3.44)  by  M  ,  and  solving  for  Mjij, 
we  obtain  from  (3,43) 

(K  * ■ p  * w M  iFo4#  iIo>  •  i3-45' 

A 

Let  us  define  an  effective  stiffness  matrix,  K  ,  by 
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K  =  K  + 


BAt 


and  an  effective  load  vector  by 

P  =  ~BAt  M^0^TT  FcP  ‘ 

Then  (3.45)  takes  on  the  form 

=  P  .  (3.46) 

After  computing  u,  ,  using  an  elimination  scheme  which  takes 
advantage  of  the  symmetric,  banded,  and  sparse  structure  of  K  ,  we 
compute  jij  and  |i^  from  (3.44):  thus, 

Hi  •  ^?{f‘i'Ho“'(r6,it2Po]  (3-47> 

and 

Fl  =  Ko  +  2At{Fo  +  Fl  }*  (3.48) 

A 

Note  that  the  matrix  K  depends  on  the  current  geometry  through  the 
local-global  element  transformation  matrices.  The  displacements, 
are  used^to  update  those  local  transformations,  which  in  turn  are  used  to 
update  K  .  Also,  the  net  load  vector,  P  ,  is  updated  based  on  the 
increment  in  external  loads  and  the  new  internal  loads  at  time,  t+At  , 
or: 
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(3.49) 


P( t+At)  =  Pfixt(t+At)  -  P1nt(t+At), 
where  the  internal  resisting  load  is  given  by 

Pint'tt4t>  '  W1'  +  %  '  <3'50> 

The  external  loads,  Pext^t+At^  *  are  a  C0m':)inat10n  user-supplied 
external  loads  and  hydrodynamic  loads  which  are  themselves  dependent  on 
the  current  deformed  state  of  the  cable. 

The  hydrodynamic  loads  are  computed 
in  terms  of  the  local  components  of  the 
velocity  of  the  fluid  flowing  past  each 
element  of  the  cable.  (For  the  purposes 
of  this  discussion  we  will  denote  the 
local  coordinates  within  each  element  by 
(x,y,z)  as  shown  in  Figure  3.6).  The 
x  axis  of  the  cable  element  lies  along 
the  locus  of  shear  centers.  With 
respect  to  this  locus  the  aerodynamic 
center  (A.C.)  and  the  modulus-weighted 
centroid  (c.g.)  are  described  by  coor- 
di nates  (Jac,2ac)  and  (ycg,zcg)  . 
respectively.  For  each  cable  element 
there  is  a  transformation  matrix,  Eg  , 
which  transforms  the  local  coordinates 
to  the  global  system;  this  is  indicated 
in  Figure  3.7. 

Figure  3.6 

Local  Coordinate  for  Element  B 
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z 


e2 


If  are  unit^ 

vectors  along  the  local  x,y,z 
axes,  then  Eg  is  given  by 


Figure  3.7 

Sketch  Showing  Global/Local  Coordinates 


ell  e12  e13 
EB  =  e21  e22  e23 

e31  e32  e33 


(3.51) 


where  ei  =  (eji »e2i *e3i )T  ls  in  global  coordinates. 


Let  V  be  the  velocity  vector  of  the  fluid  flow  relative  to  the 
cable.  Let  the  velocity  have  global  coordinates,  V  =  (V,  and 

let  the  corresponding  local  coordinates  (at  element  B)  be  expressed  by 


vi  vi 

4  V2  • 


(3.52) 


Clearly,  V^e^  is  the  local  spanwise  component  of  V  ,  the  velocity; 


U  =  V  -  Vjej 


(3.53) 


is  the  local  normal  component  of  V  . 


The  local  angle  of  attack,  a  ,  is  given  by 


-1  ^3 

a  =  tan  1  —  . 


(3.54) 


•°*r-  -  n'ecwmii 


The  local  values  of  hydrodynamic  lift  (per  unit  length),  normal  drag 
(per  unit  length),  tangential  drag  (per  unit  length),  and  twisting  moment 
(per  unit  length)  about  the  hydrodynamic  center  are  given,  respectively, 
by 

L  =  (CL  -a)Q’C 
a 

dn  ■  V°'C  • 

DT  =  CDt"5*-C  • 
and 


MHC  ■  C«hc-Q’c2  ’  (3'55) 

where 

C  =  section  chord  length. 


«  section  lift-curve  slope, 
a 

Cn  =  normal  drag  coefficient, 

UN 

Cq  =  tangential  drag  coefficient, 
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and 


C,.  =  moment  coefficient  about  the 

MHC 

hydrodynamic  center. 


The  local  components  (F^.F^.F^)  of  these  hydrodynamic  forces  are 


F1  =  sin  (V1)Dt  , 


F^  =  COS  a  -  L  Sin  a  , 


and 


F_  =  D„  sin  a  +  L  cos  a  . 
3  N 


(3.56) 


Consequently  the  global  components  (Fj.F2.F2)  are  given  by 


V 

.  A  _ 

F1 

* 

F2 

=  eb 

F2 

A 

.F3_ 

.  F3_ 

(3.57) 


Since  the  hydrodynamic  loads  act  at  the  hydrodynamic  center  of  the 
section,  they  contribute  to  the  moment  about  the  “x  axis;  consequently. 


Mx  =  (yHCF3‘zHC^  +  MHC  * 


(3.58) 


The  global  components  (M  ,  M  ,  M  )  of  this  moment  are  given  by 

x  y  z 
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r  A  _ 

M- 
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M 

=  ED 

0 

y 

B 

M 

0 

z  J 

(3.59) 


Of  considerable  interest,  for  an  interpretation  of  the  cable's 
behavior,  are  the  trail  angle,  $  ,  and  kite  angle,  8  .  If  we  regard  the 
global  x-z  plane  as  the  plane  of  the  tow,  then  we  define  <f>  and  8 
from 


cos  <j>  =  en  , 
sin  sin  e  =  e21  , 


and 


sin  <J>  cos  0  =  e31  .  (3.60) 

2  2 

We  note  that  for  e3^+  =0  the  kite  angle  is  not  defined,  jcosa|=l  , 

thus  the  cable  is  horizontal. 

We  remark  that  with  this  formulation  we  can  easily  take  into  account 
the  pressure  distribution-induced  distortions  of  those  cable  element 
cross-sections  which  experience  a  local  flow  at  a  non-zero  angle  of 
attack.  That  is,  using  the  results  of  the  two-dimensional  analysis 
(c.f.  section  2  above)  we  let  the  hydrodynamic  coefficients  depend  on  the 
local  flow  velocity  and  the  instantaneous  angle  of  attack.  Thus,  e.g., 

Cn  h  cn  (<*,V)  ;  (3.61) 

UN  UN 


-60- 


such  a  function  could  be  given  in  tabular  form  or  introduced  through  the 
use  of  interpolating  polynomials.  This  feature  is  not  incorporated  into 
the  present  version  of  the  three-dimensional  program,  TOWLINE,  due  to 
constraints  on  the  computer  budget  for  the  present  study.  An  extensive 
series  of  two-dimensional  HYDROSAP  analyses  would  be  required  in  order  to 
construct  the  tables  for  each  individual  cable  cross-section,  under 
investigation. 

Notwithstanding,  some  samples  of  the  TOWLINE  computations  are  dis¬ 
cussed  in  section  3.3.2  below. 


3.3  Sample  Computations 

3.3.1  Cross-Section  Mechanical  Properties 

We  will  seek,  for  this  example,  the  cross-sectional  mechanical 
properties  of  an  NACA  0020  cable  section,  having  a  three-inch  chord.  The 
section  has  a  glass  fibre  strength  member  (core)  in  the  nose,  and  a  flex¬ 
ible  afterbody  made  of  rubber  (tail),  as  shown  in  Figure  3.8.  Let  us 
suppose,  in  addition,  that  the  core  extends  aft  of  the  section's  nose  to 
the  20%  chord  station.  The  core  material  is  isotropic,  with  a  Young’s 
modulus  of  E=9.4xl0^  psi;  also,  suppose  it  has  a  shear  modulus 

5 

G=3. 615x10  psi.  The  material  constants  for  the  tail  are  assumed  to  be 
E=3xl03  psi,  and  3G=1.154xl03  psi. 

We  triangulate  the  cross-section  as  shown  in  Figure  3.8.  Fifty- 
eight  triangular  elements,  defined  by  46  node  points,  make  up  the  finite 
element  mesh.  We  expect,  due  to  the  large  ratios  of  the  material  prop¬ 
erties  (for  the  nose  compared  to  those  of  the  tail)  that  the  mechanical 
properties  of  the  section  will  essentially  be  those  of  the  nose.  Our 
expectations  are  confirmed  by  the  results  shown  in  Table  3.2,  where  we 
list  the  mechanical  properties  of  the  heterogeneous  section,  the  core 
alone,  and  for  a  homogeneous  section  having  the  material  properties  of 
the  core  throughout  the  entire  section.  Clearly,  the  concentration  of 


Table  3.2 

Cross-sectional  Properties  of  the  NACA  0020 
Cable  Fairing 


MECHANICAL 

PARAMETER 

THE 

HETEROGENEOUS 

SECTION 

— 

THE 

HOMOGENEOUS 

SECTION 

Modulus-weighted 
Centroid  (in). 

y 

z 

0.35136 

0. 

0.34996 

0. 

1.2623 

0. 

Shear  Center  (in.) 

ys 

zs 

0.36018 

0. 

0.35898 

0. 

1.1129  . 

0. 

Modul us-weighted 
Bending  Stiffness 

(in4) 

ETyz/E*  (1) 

nzz/E* 

4.745xl0~3 

6.925xl0"3 

4.738xl0'3 

6.409xl0~3 

2.521xl0“2 

6.044xl0_1 

Modul us-weighted 
Torsional  Rigidity 

(in4) 

gj/g*  (2) 

8.978xl0"3 

8.938xl0-3 

1.085X10'1 

Modul  us-weighted 

2 

Area  (in  ) 

ea/e* 

0.25504 

0.25470 

1.2252 

(‘)  E*  =  9.4xl06  psi  (2)  G*  =  3.62xl06  psi 
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Young's  modulus  and  shear  modulus  in  the  core  material  enables  one  to 
essentially  ignore  the  contributions  to  stiffness,  etc.,  due  to  the  tail 
material . 

3.3.2  The  Towline's  Response 

The  TOWLINE  computer  program  is  designed  to  determine  the  (static) 
equilibrium  position  of  an  integrated  (faired)  towline  under  steady  tow¬ 
ing  conditions.  This  is  to  be  accomplished  by  simulating  the  dynamic 
response  of  the  tow  cable  to  the  hydrodynamic  forces  induced  by  the  tow¬ 
ing  platform  and  the  cable's  three-dimensional  dynamic  response.  The 
TOWLINE  code  is  an  adaptation  of  a  program  written  by  I ITRI  [18]  to  study 
the  large  displacements  and  deformations  incurred  in  vehicle  collisions. 

In  this  section  we  present  the  results  of  computations  on  a  model 
cable  problem  which  was  used  for  program  development.  In  order  to  remain 
within  the  constraints  on  computer  resources,  the  model  problem  has  only 
nine  elements  and  has  a  higher  stiffness-to-length  ratio  than  would  be 
found  in  an  operational  towline.  The  results  are  intended  only  to  demon¬ 
strate  the  qualitative  behavior  of  faired  tow  cables;  and,  to  demonstrate 
the  capabilities  of  the  TOWLINE  program. 

Suppose  that  a  108-inch  faired  cable  is  suspended,  vertically,  in  a 
264  in/sec  steady  fluid  flow,  and  that  the  cable  is  subject  to  a  2.64  in/ 
sec  cross-current.  By  these  conditions,  the  cable's  cross-sections  have 
an  initial  angle  of  attack  of  0.575°  with  respect  to  the  remote  resultant 
flow.  The  cable's  upper  end  (point)  is  "fixed"  in  space,  but  free  to 
swivel.  No  depressor  is  fixed  to  the  cable's  lower  endpoint,  which  is 
completely  unconstrained  in  space.  The  initial  global  configuration  of 
the  cable  is  shown  in  Figure  3.9.  Although  the  cable  is  of  a  composite 
construction,  we  assume  that  effective  homogeneous  mechanical  properties 
have  been  determined  (using  program  SHCENT  as  discussed  in  section  3.1, 
above).  These  mechanical  properties  of  the  cable  are  summarized  as  fol¬ 
lows: 
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Local  node 2 
Global  node  10 


-4  2  4 

cable  density;  3.7x10  Ibf-sec  /in 


_2 

•  weight  in  water;  8.46x10  Ibf/in 

—  5 

•  modulus  weighted  area;  EA  =  8x10  Ibf 

•  position  of  modulus  weighted  centroid;  0.5  in.  aft  of  leading 
edge 

•  modulus  weighted  moments  of  inertia: 


EIyy  =  lxlO4  Ibf- in2  , 

IT  2z  =  4.67xl05  Ibf- in2  , 
EIyz  =  0  (a  symmetric  section) 


and 


•  torsional  rigidity; 


GJ  =  1.64xl04  Ibf- in2  . 

We  assume  that  (with  negligible  error)  that  the  shear  center  and  modulus- 
weighted  centroid  coincide.  The  hydrodynamic  properties  of  the  cable  are 
parameterized  by: 


Chord  length  of  the  cross-section;  c=2  in. 


We  model  the  cable  using  nine  beam  elements,  each  of  which  is  12 
inches  long.  A  time  step  of  At=.01  sec.  has  been  used  for  the  integra¬ 
tion  procedure.  In  Figures  3.10A  and  3.10B  we  show  the  successive  states 
of  deformation,  for  the  cable,  during  the  first  ten  integration  steps. 

We  remark  that  the  "equilibrium"  position  of  the  cable  has  (of  course) 
not  been  attained  in  only  ten  time  steps. 

Referring  to  Figure  3.10A,  we  note  that  under  the  influence  of  the 
264  in/sec.  tow  velocity,  the  cable  begins  to  move  aft;  and,  the  largest 
trail  angles  are  near  the  tow  point,  as  expected.  We  see,  however,  that 
between  the  eighth  and  tenth  integration  steps  a  wave-like  motion  has 
developed.4  While  the  free  end  of  the  cable  continues  moving  aft,  cable 
elements  9,  8,  7,  (and  possibly  6)  (nodes  1,  2,  3,  4,  5))  are  moving  for¬ 
ward.  Superimposed  on  this  swinging  motion,  aft,  is  a  higher  frequency 
vibratory  motion  in  the  vertical  direction.  This  is  typified  by  the  x-z 
position  of  node  10. 

From  Figure  3.10B  we  see  that  the  cable,  under  influence  of  the 
-2.64  in/sec  transverse  current,  is  "fluttering"  about  the  x-z  plane  of 
the  tow.  In  integration  steps  1  through  4  the  cable  moves  --  nearly  as  a 
rigid  body  —  with  the  cross-current.  Inertia  and  angle  of  attack 

4  This  "whipping"  effect  was  even  more  pronounced  in  longer  cables;  in 
many  instances  dominating  the  responses  and  sending  the  program  to 
(numerical)  divergence. 
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effects  then  swing  the  cable  back,  in  opposition  to  the  cross-current 
(steps  6  through  8).  By  integration  step  10  the  cycle  is  nearly  com¬ 
plete,  with  the  cable  approaching  the  plane  of  the  tow  and  moving  with 
the  cross-current.  Note  that  the  scales  used  to  plot  the  node  point 

.  p 

displacements  vary:  inches  times  10  for  the  x-  and  y-displacements, 

_3 

and  "inches  x  10  "  for  the  z-displacements  (vertical). 

In  Figure  3.11  we  plot  the  time  history  of  the  average  angle  of 
attack  along  the  cable.  That  the  average  angle  of  attack  is  a  reasonable 
indicator  of  local  behavior  is  shown  by  Figure  3.12.  The  twist  of  the 
cable  is  only  marginally  less,  near  the  cable's  attachment  point,  than  at 
the  free  end.  The  cross-current  of  2.64  in/sec.  corresponds  to  an  ini¬ 
tial  angle  of  attack  of  0.573°  with  respect  to  the  resultant  flow. 
Twisting,  induced  by  the  aft  location  of  the  hydrodynamic  center,  with 
respect  to  the  shear  center,  causes  the  cable's  angle  of  attack  to  vary 
from  the  initial  (0.573°)  to  -1.45°  (at  step  6),  then  back  to  2.16°  at 
step  8,  and  to  4.72°  at  step  10.  This  oscillation  in  the  torsional  dis¬ 
placement  of  the  cable  contributes  to  the  oscillatory  kiting  (transverse 
motion).  Moreover,  even  from  the  limited  duration  of  these  sample  compu¬ 
tations,  it  is  clear  that  both  the  magnitude  of  the  kiting  displacement, 
and  the  angle  of  attack,  are  growing,  even  though  the  hydrodynamic 
center's  position  —  aft  of  the  shear  center  —  gives  a  "statically 
stable"  configuration. 

We  note  once  again  that  the  physical  parameters  of  the  cable,  anal¬ 
yzed  above,  are  not  representative  of  operational  cable  systems.  This 
sample  cable  has  a  very  high  stiffness-to-length  ratio,  compared  to  an 
operational  system.  Accordingly,  the  frequencies  and  the  magnitudes  of 
the  oscillations,  exhibited,  are  not  representative  of  operational 
cables.  The  complete  analysis  of  an  operational  systen  could  not  be  per¬ 
formed  during  the  present  developmental  effort  due  to  limitations  on  the 
computer  budget. 
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4.0  CONCLUSIONS  AND  RECOMMENDATIONS 


The  initial  phases  in  the  development  of  a  computational  program,  to 
be  used  for  studying  faired,  flexible  towing  cables,  has  been  completed. 
In  the  first  phase  of  this  task  effort,  the  chordwise  deformations,  pro¬ 
duced  on  a  flexible  fairing  profile,  in  a  steady,  two-dimensional  viscous 
flow,  were  obtained  using  the  software  module  HYDROSAP.  This  is  now  an 
operational  system  which  can  be  used  on  a  production  basis  to  describe 
the  deformations  of  flexible  fairings  under  hydrodynamic  loadings. 

In  this  document  the  mathematical  and  physical  foundations  upon 
which  the  HYDROSAP  software  is  based  were  summarized.  In  addition,  the 
capabilities  of  the  system  were  illustrated  through  the  presentation  of 
typical  results  acquired  from  its  application  to  several  model  cable  pro¬ 
files.  Even  though  it  was  not  emphasized  here,  it  should  be  noted  that 
this  program  is  a  fully  developed  tool  for  analysis;  developed  to  the 
extent  that  it  includes  highly  automated  input  generator  schemes  which 
make  it  particularly  attractive  to  the  analyst-user.  In  this  regard  the 
user  is  relieved  of  the  troublesome  task  of  building  lengthy  and  compli¬ 
cated  input  files.  Only  a  minimum  of  control  and  descriptive  input 
information  needs  to  be  provided  in  order  to  run  this  program. 

During  the  present  task  significant  steps  were  made  toward  the  com¬ 
pletion  of  a  second  phase  in  this  program's  development  effort.  Exten¬ 
sions  to  the  two-dimensional  computational  capabilities,  moving  toward  a 
full  three-dimensional  capability,  have  been  embodied  —  in  preliminary 
form  --  into  two  additional  computer  programs,  SHCENT  and  TOWLINE. 

As  we  noted  in  section  3.0,  the  three-dimensional,  dynamic  behavior 
of  a  faired  cable  system,  which  is  modelled  in  TOWLINE,  uses  a  small 
strain,  large  displacement  incremental  finite  element  algorithm.  For 
this  model,  the  traditional  twelve  degree  of  freedom  beam  element,  of 
linear  structural  analysis,  has  been  adapted  to  the  present  incremental 
formulation.  Of  course,  the  use  of  this  beam  element,  to  discretize  the 
three-dimensional  equations  of  motion  for  the  cable,  requires  a  knowledge 
of  certain  of  the  mechanical  properties  (bending  and  torsional  rigid- 
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ities;  shear  center  location)  for  the  cable's  cross-sections.  Neces¬ 
sarily,  the  mechanical  properties  of  heterogeneous,  streamlined  cross- 
section  profiles,  used  to  describe  faired  cables,  must  themselves  be  com¬ 
puted.  We  have  discussed  in  this  docunent  how  the  SHCFNT  program  per¬ 
forms  that  function.  Although  it  is  not  yet  as  "user-friendly"  as  the 
HYDROSAP  module,  SHCENT,  too,  is  available  for  production  use. 

Even  though  the  TOWLINE  program,  itself,  is  operational  —  as  we 
have  discussed  in  section  3.0  above  —  it  is  not  yet  ready  for  computa¬ 
tional  production  ,.ork.  First,  its  input  requirements  are  rather  lengthy 
and  cunbersome.  No  effort  has  been  expended  in  giving  it  an  automated 
"user-friendly"  input  capability.  Second,  we  have  found,  through  an 
extensive  series  of  check  runs,  several  undesirable  performance  charac¬ 
teristics  which  would  detract  from  its  use  in  the  eyes  of  a  routine  pro¬ 
duction  user. 

Due  to  the  combination  of  an  extreme  slenderness  ratio  (span/chord), 
as  is  typical  of  operational  towline  systems,  and  the  use  of  an  incre- 
mental  dynamic  solution  technique,  in  following  the  towline  to  a  desired 
equilibria  state,  we  have  found  that  the  TDWLINE  algorithm  demands 
inordinately  small  time  steps.  These  small  time  steps  are  necessary  to 
follow  the  complicated  high-frequency  dynamic  responses  which  are  appar¬ 
ently  exhibited  by  these  cable  systems.  High  frequency  oscillations,  in 
addition  to  cable  "whipping",  were  observed  in  several  of  the  sample  com¬ 
putational  runs.  Since  the  required  integration  time  step  is  of  the 
order  of  micro- seconds,  it  is  prohibitively  expensive  to  follow  a  motion 
for  a  physically  meaningful  duration  (for  several  minutes,  say).  Realis¬ 
tically,  the  time  propagation  of  numerical  errors,  inherent  to  any  compu¬ 
tational  algorithm,  would  be  of  serious  concern  for  any  run  of  such  mag¬ 
nitude. 

Based  on  our  experience  with  the  current  version  of  TOWLINF,  then, 
we  recommend  that  serious  attention  be  given  to  the  development  of  a 
non-linear  cable  finite  element.  Any  such  new  element  formulation  would 
be  tailored  to  the  modelling  of  the  structural  response,  for  structural 
members  —  such  as  cables  --  with  very  high  slenderness  ratios.  However, 
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unlike  existing  cable  elements  which,  at  most,  are  capable  of  undergoing 
axial  extensions,  our  proposed  cable  element  would  also  resist  both  bend¬ 
ing  and  torsion.  This  new  cable  elenent  formulation  should,  of  course, 
not  be  based  on  simple  models  from  the  linearized  theory  of  strength  of 
materials,  but  must,  instead,  incorporate  the  more  recent  advances  in 
thin  rod  theories  [19,  20,  21]. 

Therefore,  a  general  recommendation,  based  on  results  acquired  from 
this  work,  is  that  the  presently  used  (cable  finite  element)  model  should 
be  replaced  by  a  more  representative  one.  With  such  a  properly  designed 
element  it  is  anticipated  that  many  of  the  present  systan’ $  computational 
difficulties  will  be  overcome.  It  is  expected  that  when  this  new  element 
model  is  introduced  into  (say)  the  TOWLINE  program,  it  will  provide  an 
efficient  and  useful  tool  for  engineering  design  and  analysis.  Concur¬ 
rent  with  the  development  of  this  new  finite  elenent,  it  is  also  recom¬ 
mended  that  the  present  computational  efforts  be  continued  to  their 
intended  conclusion.  Then,  the  completed  system  can  and  should  be  placed 
into  operation  as  a  design  and  study  tool  for  high-speed,  flexible  tow- 
line  systems. 
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APPENDIX  A  Hydrosap  Computer  Program 
A.l  Program  Description 

The  Hydrosap  computer  program  computes  the  deformed  configuration  of 
a  two-dimensional  "flexible"  fairing  at  angle  of  attack  in  a  steady  vis¬ 
cous  fluid  flow.  The  program  consists  of  three  segments:  HYDRO,  for  the 
"fluids"  computations;  NONSAP,  for  the  mechanical  response  computations, 
and  DRIVER,  an  executive  sub-program  which  initializes  and  controls  the 
computations.  The  NONSAP  segment  consists  of  routines  selected  from  the 
University  of  California  Non-linear  Structural  Analysis  Program  while  the 
HYDRO  segment  is  an  adaptation  of  the  NASA  Multi-component  Two- 
dimensional  Viscous  Airfoil  Program.  Since  we  have  discussed  the  details 
of  HYDRO  and  NONSAP  operations  elsewhere  [1],  and  since  the  user  "sees" 
only  the  front-end  executive,  DRIVER,  we  discuss  in  this  Appendix  only 
the  DRIVER  input  deck. 

A. 2  DRIVER  Input  Deck 

The  HYDROSAP  input  cards  are  identified  with  Tape  5  (usually  the 
system  inpu*-  stream).  The  input  stream  consists  of  six  distinct  card 
groups,  not  all  of  which  are  present  in  a  single  run.  Only  the  first  two 
card  groups  must  be  present;  the  presence  or  absence  of  the  other  four 
groups  is  determined  from  values  assigned  by  the  user  to  the  NONSAP  con¬ 
trol  variables  SHRTIN,  HYDSIM,  SAPSIM,  and  N.  All  of  the  control  vari¬ 
ables  are  defined  below.  The  input  stream  groups  are  as  follows: 


CARD  GROUP  1  Title  Card 

Card  1.1  Title  of  Run 

T 

Up  to  an  80  character  title  for  the  run. 
(This  card  must  appear.) 
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CARD  GROUP  2  Control  Variables 

Card  2.1  $DRV 

12 

Header  card  for  NAMELIST  /DRV/;  signals  the  beginning  of 
the  control  input.  (This  card  must  appear.) 

Definitions  of  the  DRV  NAMELIST  variables  are  as  follows: 

A. 2.1  Control  Variable  Definitions 


The  HYDROSAP  Control  Variables  (Card  Group  2  of  the  HYDROSAP  Input)  fall 
into  two  groups:  logical  variables  and  arithmetic  variables.  The  logical 
variables  are  true/false  switches  used  to  activate  the  various  program  options. 
The  arithmetic  variables  are  used  to  set  values  for  the  integer  and  real 
variables  which  control  the  flow  of  computations  (examples  are:  iteration 
counters,  output  units,  mesh  fineness,  and  convergence  tolerance),  or  define 
the  physical  parameters  of  the  problem  under  study  (such  as  flow  conditions 
and  material  properties).  Detail  definitions  of  the  control  variables  are 
as  follows: 


VARIABLE 

CHKPNT 

(CHeckPoiNT) 

CORE 


DEBUG 


TYPE  DEFAULT  INTERPRETATION 


LOGICAL  .FALSE.  =.TRUE.  Produce  checkpoint 

tape  (FT14F001) 

=. FALSE.  Checkpointing  dis¬ 
abled. 

LOGICAL  .FALSE.  =.TRUE.  Fairing  has  a  cable 

(core)  material 
=. FALSE.  No  distinct  cable 
(core)  present 
(SHRTIN= .TRUE,  only) 

LOGICAL  .FALSE.  =.TRUE.  Special  print  in 

LOADER:  HYDRO  and 
NONSAP  output  pro¬ 
vided  at  each  HS 
iteration  (itera¬ 
tion  =  outer  loop) 

=. FALSE.  HS  print  produced  at 
start  of  run  only. 

No  debug  output. 
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VARIABLE 

TYPE 

DEFAULT 

INTERPRETATION 

ECHO 

LOGICAL 

.TRUE. 

=.TRUE.  Echo  of  HYDRO  and 

NONSAP  input  decks 
produced. 

=. FALSE.  Suppress  "echo"  of 
files. 

HYDSIM 

(HYDro  SIMulation) 

LOGICAL 

.FALSE. 

=.TRUE.  Program  operates  in 
the  HYDRO  simulation 
mode  (SAPSIM  must  be 
.FALSE.) 

=. FALSE.  Signifies  a  request 
for  the  full  HYDROSAP 
run;  or  NONSAP 
simulation,  depend¬ 
ing  on  the  value 
assigned  to  SAPSIM. 

POST 

LOGICAL 

.FALSE. 

=.TRUE.  The  Graphics  Post¬ 
processor  output  is 
directed  to  Unit  88. 

=. FALSE.  No  operation  requested. 

RESTRT 

LOGICAL 

.FALSE. 

= .TRUE .  Read  checkpoint  tape 
and  restart  at  time 
TSTART 

*. FALSE.  Restart  disabled. 

SAPSIM 

LOGICAL 

.FALSE. 

=.TRUE.  Program  is  to  operate 
in  NONSAP  simulation 
mode  (HYDSIM  must  be 
set  to  .FALSE.) 

=. FALSE.  Full  HYDROSAP  run; 

or  HYDRO  simulation, 
depending  upon  value 
assigned  to  HYDSIM. 

SHRTIN 

(SHoRT  form  INput) 

LOGICAL 

.TRUE. 

=.TRUE.  Signifies  the  short 
form  input.  An  auto¬ 
matic  input  file  gen¬ 
eration  occurs  (for 
foil  geometry  and 
state) . 

=. FALSE.  HYDRO  and  NONSAP  in¬ 
put  decks,  separated 
by  999999999,  must 
follow  $DRV-$END  on 
SYSIN  file. 
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VARIABLE 

TYPE 

DEFAULT 

INTERPRETATION 

SKIN 

LOGICAL 

.FALSE. 

=.TRUE.  Denotes  the  fairing 
cross-section  has  a 
skin  attached. 

=. FALSE.  No  skin  member  is 
presumed. 

TAIL 

LOGICAL 

.TRUE. 

=.TRUE.  The  fairing  has  a 
fairing  "tail"  mem¬ 
ber. 

=. FALSE.  No  fairing  tail  is 

assumed  (CORE  must  be 
set  to  .TRUE,  in  this 
case) . 

TPRINT 

(Terminal  PRINT) 

LOGICAL 

.TRUE. 

=.TRUE.  Terminal  summary  out¬ 
put  (error  msgs.  and 
job  flow) 

=. FALSE.  No  terminal  print. 

INTORD 

(INTerpolation  ORDer) 

INTEGER 

3 

Denotes  the  interpolation 
order  for  the  pressure  dis¬ 
tribution  conversion  from  the 
HYDRO  surface  to  the  NONSAP 
surface. 

INTORX 

(INTerpolation  ORder  X) 

INTEGER 

1 

Interpolation  order  to  be 
assigned  for  the  displacement 
conversion  from  the  NONSAP 
surface  to  the  HYDRO  surface. 

I  PANEL 

INTEGER 

2 

I PANEL  is  equivalent  to  HYDRO 
control  variable  IOP. 

IPAPER 

INTEGER 

6 

Unit  number  assigned  to  the 
HYDROSAP  DRIVER  module  output. 

IPRINT 

INTEGER 

=NSTEPS 

Denotes  the  NONSAP  print  fre¬ 
quency  within  HYDROSAP 
ITERATION  #1. 

ISMTH 

INTEGER 

0 

ISMTH  is  equivalent  to  HYDRO 
control  variable  ISMO. 

ITRSWT 

(ITeRation  SWiTch) 

INTEGER 

8 

The  number  of  HYDROSAP  itera¬ 
tions  using  HYDRO  potential 
flow  solution  mode  only. 

IWRT0 

INTEGER 

6 

Denotes  the  NONSAP  output 

UNIT. 
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VARIABLE 


TYPE 


DEFAULT 


INTERPRETATION 


IWRT1 

INTEGER 

6 

Denotes  the  HYDRO  output  unit. 

LOOPMX 

(LOOPMaXimum) 

INTEGER 

1 

Signifies  the  maximum  number 
of  HYDROSAP  iterations  (outer 
loops)  allowed. 

MESHC 

(MESHChordwise) 

INTEGER 

1 

Denotes  the  number  of  chord- 
wise  interior  stations  used 
to  generate  the  finite  ele¬ 
ment  mesh,  when  SHRTIN  is 
set  to  .TRUE. 

MESHT 

(MESHThickness-wise) 

INTEGER 

1 

Denotes  the  number  of  lateral 
interior  stations  used  to  aen- 
erate  the  finite  element  mesh 
when  SHRTIN  is  set  to  .TRUE. 

N 

INTEGER 

0 

N>0  Signifies  the  number  of 

upper  surface  coordinate 
pairs  to  be  input  through 
NAMELIST  $FAIRNG-$END,  fol¬ 
lowing  NAMELIST  $DRV-$END, 
on  SYSIN. 

N=0  Means  no  fairing  profile 
is  input. 

N<0  A  pointer  for  profile 
library  routine  |N j . 

Thus,  N= -2  reads  profile 
points  from  subroutine 
FOIL02. 

(Applicable  when  SHRTIN  set  to 
.TRUE.,  only.) 

NSTEPS  INTEGER  20  Denotes  the  number  of  (quasi  - 

static)  load  steps  to  be  used 
in  applying  the  field  pressure 
load  to  a  fairing,  in  NONSAP, 
during  HYDROSAP  iteration  #1. 

Note:  NSTEPS  should  be  increased  when  the  message  "OUT-OF-BALANCE 
LOADS  EXCEED  INCREMENTAL  LOADS"  is  obtained. 


ALPHA 

REAL*8 

0.D0 

Angle  of  attack  (deg. ) 

ECORE 

REAL*8 

1.D06 

Young's  modulus  for  the  core 

material  (psi). 

ESKIN 

REAL*8 

1  .D04 

Young's  modulus  for  the  skin 

material  (psi). 
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VARIABLE  TYPE 

DEFAULT 

INTERPRETATION 

ETAIL 

REAL*8 

1  .D04 

Young's  modulus  for  the  fair¬ 
ing  tail  material  (psi). 

PCORE  1 

1 

.3 

Poisson  ratio  for  the  core. 

PSKIN 

REAL*8 

0 

skin,  and  tail , 

PTAIL  j 

r 

.3 

respectively. 

PZERO 

REAL*3 

I4.7D0 

Static  pressure  loading  (psi). 

RHO 

REAL*8 

2.D0 

H20  density  (slugs/ft3). 

TOL 

REAL *8 

lxl0-2 

Denotes  the  HYDROSAP  conver¬ 
gence  tolerance  level.  Con¬ 
vergence  is  achieved  when  rms 
displacement  increment 
<(rms  displacement  )  *  TOL. 

TSKIN 

REAL*8 

.005 

Skin  thickness  (inches) 

(applies  only  if  SKIN  is  set 
to  .TRUE.). 

TSTART 

LOGICAL 

0.D0 

Applies  only  when  RESTRT= 

(Time  of  reSTART) 

.TRUE.  HYDROSAP  will  restart 
at  checkpoint  for  time  TSTART. 

VFRSTR 

REAL*8 

22.  D0 

Freestream  velocity  (ft/sec) 

VFRSTR 

REAL *8 

22.  D0 

Free  stream  velocity  (fps). 

Y1CORE 

REAL*8 

0.D0 

Chordwise  location  of  the  core 
starting  point. 

Y2CORE 

REAL*8 

0.D0 

Chordwise  location  of  the  core 
ending  point. 
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Card  2.2 


SHRTIN=T,  •  •  • 

123 

Input  as  many  cards  as  required  to  override  default  control 
variable  values. 

Card  2.n 

Card  2.n+l  $END 
12 

End  card  for  NAMELIST  /DRV/;  signals  the  termination  of 
the  namelist  string.  (This  card  must  appear.) 


CARD  GROUP  3  Cable  Fairing  Profile  Namelist  (optional) 

Card  Group  3  is  present  only  if  SHRTIN=T  and  N>0  in 

NAMELIST  /DRV/. 

Card  3.1  $FAIRNG 
12 

Header  card  for  NAMELIST  /FAIRNG/ . 

Card  3.2  XX(1 )  =  • • • ,XX(N)=- • • , 

123 


Card  3.m  ZZ(1 )=• • • ,ZZ(N)=* • • , 

123 

{ (XX ( I )  ,ZZ( I ) ) »  1=1 , - • * N >  are  the  profile  coordinates,  in 

inches,  for  the  upper  surface  of  the  symmetric  cable  fair 

ing  to  be  modeled.  N  is  set  in  NAMELIST  /DRV/. 

Card  5.m+l  $END 
12 

End  card  for  NAMELIST  /FAIRNG/. 


CARD  GROUP  4  HYDRO  Input  Deck  (optional) 

Card  Group  4  is  present  only  if  SHRTIN=F  and  SAPSIM=F  in 
NAMELIST  /DRV/. 

Card  4.1  First  card  of  HYDRO  input  deck 
•  }  see  section  3.3 
Last  card  of  HYDRO  input  deck. 


Card  4.p 


CARD  GROUP  5  Spacer 

Card  5.1  99999999 

12345678 

A  spacer  card  to  be  inserted  between  HYDRO  and  NONSAP  in¬ 
put  decks.  (This  card  must  be  present  whenever  SHRTIN=F.) 


CARD  GROUP  6 


Card  6.1 

Card  6.q 


NONSAP  Input  Deck  (optional) 

Card  Group  6  is  present  only  if  SHRTIN=F  and  HYDSIM=F  in 
NAMELIST  /DRV/. 

First  card  of  the  NONSAP  input  deck. 

•  }  see  section  4.3 
Last  card  of  the  NONSAP  input  deck. 
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APPENDIX  B.  SHCENT  Computer  Program 
B.l  Program  Description 

The  SHear  CENTer  (SHCENT)  program  computes  the  St.  Venant  warping 
function  for  an  arbitrarily  shaped,  heterogeneous,  orthotropic  cross- 
section.  The  SHCENT  program  was  adapted  from  a  code  developed  by  E. 
Baumgartner  [15].  Using  the  warping  function,  the  program  computes  the 
section's  torsional  rigidity  and  the  coordinates  of  its  shear  center. 
Incidental  to  the  warping  function  computations,  the  program  also  com¬ 
putes  the  bending  rigidities  of  the  section,  the  modulus  weighted  area 
for  the  section,  and  the  location  of  the  modulus  weighted  centroid. 
Written  in  FORTRAN,  the  program  accepts,  as  input,  the  topology  of  fi 
(the  section  shape)  and  its  triangulation,  as  well  as  the  material  con¬ 
stants  for  each  element  of  the  triangulation. 

B.2  Problem  Size  Limitations 


The  current  double  precision  version  of  the  program  runs  in  224k 
bytes  of  core  on  an  IBM  system  360/91.  With  this  program  size,  the  tri¬ 
angulation  of  52  is  limited  to: 

•  maximum  number  of  triangular  elements  -  140 

•  maximum  number  of  nodal  points  -  150 

•  maximum  number  of  different  materials  comprising  the 
cross-section  -  20 

•  maximum  semi-bandwidth  for  the  global  stiffness  matrix  -  100. 

All  significant  variables  are  passed  among  the  various  subroutines  of 
SHCENT  via  blank  common.  The  problem  size  limitations  (above)  are  trivi¬ 
ally  eased  by  altering  the  dimensions  of  the  variables  in  blank  common. 
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B.3  SHCENT  Input  Deck 


Program  SHCENT  reads  input  from  FORTRAN  Unit  5  using  list-directed 
(free-field)  format.  The  programmer  may  comment  his  input  deck  by  plac¬ 
ing  an  asterisk  (*)  in  the  first  column  of  the  comment  record.  A  rudi¬ 
mentary  mesh  generating  option  is  available  for  automatically  computing 
model  coordinates  and  element  definitions  for  "regular"  meshes.  For 
this,  element  nodes  must  be  "numbered"  counterclockwise. 

A  description  of  the  input  deck  follows: 

Title  Card 

Card  1:  TITLE 

TITLE  -  Any  80  characters  alphanumeric  to  be  used  as  the 
title  for  a  run 

Control  Card 

Card  2:  NODNO,  NGMAT,  NGRELM,  NGRNOD,  NECHO 

NODNO  -  Total  nunber  of  nodes 

NGMAT  -  Number  of  material  property  vectors  to  be  read 

['’ELM  -  Number  of  element  definition  groups  to  be  read 

GRNOD  -  Number  of  node  point  definition  groups  to  be  read 

NECHO  -  Input  echo  flag 

=  0,  No  echo  of  input 
=  1,  Echo  input 

Material  Property  Cards 


Cards  3:  K,  (GMAT(k.j) ,j=I ,4)  (input  NGMAT  cards) 
GMAT(k,l)  -  G^z  ,  xz  shear  constant 
GMAT(k,2)  -  Gyz  ,  yz  shear  constant 


GMAT(k,3)  -  Ez  ,  Young's  modulus  orthogonal  to  the 
cross-section 

GMAT{k,4)  -  a  ,  angle  (degrees)  between  material  axes 
and  input  coordinate  system 

Element  Definition  Cards 


Cards  4:  NGR,  IL,  ILA,  IS,  IA,  JS,  JA,  KS,  KA,  NMAT 

NGR  -  Number  of  elements  to  be  defined  by  this  card 

IL  -  Element  number  of  first  element  in  the  group 

ILA  -  Increment  by  which  IL  is  to  be  increased 

(decreased  for  ILACO)  for  each  additional  element 
defined  in  this  group 

IS  -  Global  node  number  of  element  (local  node)  number 
1 

IA  -  Node  number  increment  for  IS  (cf.  ILA) 

JS  -  Global  node  number  of  element  (local  node)  number 

2 

JA  -  Node  number  increment  for  JS 

KS  -  Global  node  number  of  element  (local  node)  number 

3 

KA  -  Node  number  increment  for  KS 

NMAT  -  Material  property  vector  to  be  associated  with 
this  element  (cf.  K  of  CARDS  3) 

Node  Point  Definition  Cards 


CARDS  5:  NGR,  J,  JA,  Y(J),  DY,  Z(J),  DZ 

NGR  -  Number  of  nodes  to  be  defined  by  this  card 
J  -  Node  number  of  first  node  in  this  group 

JA  -  Increment  by  which  J  is  to  be  increased  (or 

decreased,  for  JA<0)  for  each  additional  node 
defined  in  this  group 
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Y(J)Z(J)-local  or  global  coordinates  of  node  J 
DY.DZ  -  coordinate  increments  by  which  Y(J)  and  Z(J), 

respectively  are  to  be  increased  for  each  addi¬ 
tional  node  defined  in  this  group 


End  of  Data  Card 


Card  6:  -999  ,  in  columns  1-4. 


! 

r 

* 


L 
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APPENDIX  C  TOWLINE  Computer  program 
C.l  Program  Description 

The  TOWLINE  program  computes  the  three-dimensional  dynamic  (or 
quasi-static)  response  of  a  faired  towing  cable  under  the  assumption  that 
the  cable,  while  undergoing  large  displacements,  sustains  only  small 
strains.  A  standard  three-dimensional  beam  element  formulation  is  used 
to  discretize  the  cable. 

Written  in  FORTRAN,  the  program  accepts  as  input  an  initial  config¬ 
uration  of  the  cable  and  the  mechanical  and  hydrodynamic  properties  of 
the  elements  into  which  the  cable  has  been  subdivided.  Up  to  ten  differ¬ 
ent  beam  cross-section  types  can  be  accommodated.  The  user  may  specify 
time-dependent  external  forcing  functions,  in  tabular  format,  or  supply  a 
user-written  subroutine  (FORCE)  to  define  the  forcing  functions  explic¬ 
itly.  A  checkpoint-restart  capability  makes  possible  the  resumption  of 
computations  from  a  previous  run. 

C.2  Problem  Size  Limitations 


The  current  double-precision  version  of  the  program  runs  in  498K 
bytes  of  core  on  an  IBM  System  360/91  system.  With  this  program  size  the 
discretization  of  the  cable  is  limited  to  100  nodes  (99  elements). 

C.3  TOWLINE  Input  Deck 


Program  TOWLINE  accepts  input  from  FORTRAN  Unit  5  using  formatted 
reads.  A  description  of  the  input  variables  and  card  formats  follows: 

CARD  1:  KEY,  TITLE  (19)  FORMAT  (Al,  19A4) 

KEY  -  decision  flag  for  integration  procedure 
=  E,  explicit 
=  I,  implicit 

TITLE  -  76  character  titles  for  the  run 
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CARD  2:  NNODE,  NELE,  NUMMAT, 

NUMDIS,  MXSTEP,  NDGREE, 
DELT,  NCOORD,  NSECT,  NFUN, 
IAERO 


FORMAT  (615,  EIO.6,415) 


NNODE  -  number  of  structural  nodes  (max  100) 

NELE  -  nunber  of  elements  (max  100) 

NUMMAT  -  number  of  different  material  models 
NUMDIS  -  ntmber  of  nodes  at  which  displacement,  velocity, 
or  force  are  to  be  specified,  as  functions  of 
time 

MXSTEP  -  maximum  number  of  time  steps 
NDGLEE  -  number  of  degrees  of  freedom  per  node  (=6  for 
current  version) 

DELT  -  time  step  size  (secs) 

NCOORD  -  number  of  beam  cross-section  pointer  nodes 
NSECT  -  number  of  different  beam  cross-section  types 
NFUN  -  number  of  different  time-dependent  functions  (max 

*  9) 

IAERO  -  decision  flag  for  hydrodynamic  loads 
=  0,  no  hydrodynamic  loads 
s  1,  include  hydrodynamic  loads 


MXSTEP 

NDGLEE 


NCOORD  - 
NSECT  - 


IAERO  - 


CARD  3:  RH0,V(3) 


-  fluid  density 

-  fluid  velocity  vector,  with  respect  to  cable,  in 
global  (x,y,z)  coordinates  (x-z  plane  repre¬ 
sents  plane  of  tow;  z  axis  along  gravity 
vector) 


OMIT  CARD  3  IF  IAELO  =  0. 


.*r ijmqji 


CARO  4:  KONTRL  (10) 


FORMAT  (1015) 


KONTRL  (1) 
(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

(8) 


(9) 

(10) 


CARD  5.1:  MTYPE 


-  stiffness  matrix  update  frequency 

-  not  used 

-  not  used 

-  not  used 

-  6x1000;  Newmark  -  beta 

-  not  used 

-  not  used 

-  checkpoint/restart  control 
=  0,  no  operation 

=  n,  checkpoint  every  n  steps  on  unit  25 
=  -n,  restart  and  then  make  checkpoint  every 
steps 

-  restart  checkpoint  number 

-  time  history  frequency 
<  0,  no  operation 

=  n,  make  time  history  of  motion  on  unit  22 
every  n  steps. 

FORMAT  (15) 


MTYPE  -  material  model  number 


CARD  5.2  &  5.3  E(16)  FORMAT  (8F10.0) 

E(l)  -  material  density 

E(2)  -  modulus  of  elasticity 

E(3)-E(6)-  not  used 
E(7)  -  Poisson's  ratio 

E(8)-E(16)-not  used 

NOTE:  REPEAT  CARD  GROUP  5  NUMMAT  TIMES 


CARD  6.1:  NSEG,  KLOS,  IYY,  IZZ,  J 


FORMAT  (2I5.3F10.0) 


NOTE: 


Figure 


NSEG  -  lumber  of  plate  sections  in  beam  cross-sections 
(set  NSEG=1  if  KL0S=2) 

KLOS  -  cross-section  type 

=  0,  thin  walled  closed  section 
=  1,  thin  walled  open  section 
=  2,  full  cross-section  with  user  specified  cross- 
section 

IYY  -  user-supplied  cross-section  moments  of  inertia 
IZZ  about  local  y  and  z  axes,  respectively,  of 
the  section 

J  -  user-supplied  cross-section  torsion  constant 


CARDS  6.2:  Y( I) ,  Z( I) ,  T( I) ,  1=1 .KLOS+NSEG  FORMAT  (3F10.4) 


Y(I)  -  cross-section  coordinate  of  point  at  beginning  of 
Z(I)  plate  segment  I  (cf.  Figure  c.l  below) 

T(I)  -  thickness  of  cross-section  plate  segment  I 

(1)  If  KL0S0=2,  then  Y(3),Z(3)  are  the  coordinates  of  the 

modulus-weighted  centroid  of  the  section  with  respect  to  the 
section  coordinate  system  (origin  at  the  section's  shear 


center).  T(2)  and  T(3)  are 


(2)  Input  NSEG+KLOS  cards  in  group 


c.la  Thin  Wall  Bean  Cross  Section 
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ignored  in  this  case. 
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Figure  c. lb  Full  Beam  Cross-section 
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